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ABSTRACT 
The  objective  of  this  work  was  the  evaluation  of  pseudospectral 
methods  for  their  computational  efficiency  and  applicablity  to  problems 
in  transport  phenomena.  This  was  accomplished  by  application  of 
Chebyshev  pseudospectral  methods  to  problems  from  each  major  area  of 
transport  phenomena.  The  evaluation  of  the  pseudospectral  method  was 
based  upon  comparison  to  the  analytic  solution,  if  available,  or  finite 
difference  approximation  otherwise.  It  was  found  that  discontinuities 
in  the  solution  domain  can  result  in  serious  deviations  from  the  correct 
solution;  for  example,  the  temperature  discontinuity  in  thermal 
entrance  length  problems  led  to  the  propagation  of  error  within  the 
solution.  Digital  filtering  was  used  successfully  to  damp  out 
oscillatory  behavior  in  all  cases  studied. 
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CHAPTER  I.   INTRODUCTION 

Transport  phenomena  are  pervasive  in  chemical  engineering  and 
accurate  numerical  solutions  are  requisite  for  improvements  in  process 
design  and  analysis.  However,  the  evaluation  of  these  equations  is 
often  difficult  due  to  coupling  and  nonlinearity.  This  has  given  rise 
to  many  of  the  empirical  relations  used  by  practicing  engineers  today, 
and  is  largely  responsible  for  the  development  of  the  "unit  operations" 
approach.  While  empirical  relations  can  be  used  for  process  design  it 
is  sometimes  difficult  to  predict  exactly  what  will  happen  if  a  par- 
ticular variable  is  changed.  With  the  universal  availability  of 
computers  many,  if  not  most,  linear  systems  are  now  modeled  numerically, 
allowing  optimum  operating  characteristics  to  be  determined. 

The  most  popular  numerical  scheme  for  solving  partial  differential 
equations  is  the  finite  difference  method  which  is  easily  applied  to 
many  problems.  Finite  difference  techniques  can  be  used  in  conjunction 
with  explicit,  implicit  or  hybrid  methods  depending  upon  the  nature  of 
the  equation.  Another  popular  numerical  approach  is  the  finite  element 
method  (see  Finlayson  [1980]).  While  a  great  amount  of  numerical  work 
is  carried  out  in  all  fields  of  Chemical  Engineering  some  of  the  most 
advanced  work  has  occurred  in  the  area  of  fluid  flow,  primarily  because 
of  the  nature  of  the  Navier-Stokes  equation,  and  the  problem  of 
turbulence.   Consider  the  equations: 

""l  t    8ui   1  3p_     A 

aT  +  uj  ax-  =  paf:  +  U7T-  U> 

J       i     8x. 
and 


3u. 

where  u.  is  the  instantaneous  velocity.   This  problem  is  often  attacked 

using  Reynolds  decomposition  by  letting  u  =  U  +  u';  where  U  is  the  mean 
velocity  and  u1  is  the  fluctuation  about  the  mean.  The  resulting  equa- 
tions are  time  averaged  and  simplified  by  noting  that  the  linear  terms 
in  the  fluctuating  quantities  are  zero  for  a  statistically  stationary 
process.  The  result  of  this  procedure  is  the  Reynolds  momentum 
equation: 

3ui     „    8ui         iap       i    a    ,    9ui 

aT  +  uj  alT  "  "  p  alT  +  p  alT  ^  alT  "  p  uiu  >  (3) 

It  is  the  introduction  of  -puju1.  terms  that  constitutes  the  major  dif- 
ficulty confronting  any  analyst  wishing  to  employ  this  method  of  attack 
upon  the  problem.  An  approximation  for  the  turbulent  momentum  flux  (- 
pu'u')  is  necessary  for  calculation  of  many  flows  of  interest.  This 

term  is  usually  related  to  mean  flow  velocity  gradients  multiplied  by  a 
proportionality  function,  or  eddy  viscosity.  The  result  is  in  essence 
an  empirical  curve  fitting  approach  in  which  constants  are  adjusted 
until  the  calculated  and  experimental  results  agree.  Extrapolation  from 
these  results  can  be  dangerous.  It  has  been  noted  by  many  fluid 
dynamicists,  including  Tennekes  and  Lumley  [1972]  that  gradient 
transport  models  are  not  appropriate  for  turbent  flow  with  multiple 
length  or  velocity  scales. 

There  are  alternatives  to  the  Reynolds  equation-turbulence  modeling 
approach;   these  include  sub-grid  scale  closure  models  and  direct 


numerical  simulation  of  the  Navier-Stokes  equations.  Deardorff  [1970] 
has  been  the  principal  proponent  of  sub-grid  scale  closure — in  this 
method,  the  large  scale  motions  are  explicitly  obtained  from  the  Navier- 
Stokes  equations.  The  very  small  scale  motions  are  treated 
statistically  and  in  this  manner  the  Reynolds  number  limitations  of  a 
total  numerical  simulation  are  avoided.  Orszag  [1977],  among  others, 
speaks  in  favor  of  this  approach,  while  noting  that  the  methods 
presently  being  used  to  handle  the  dissipative  eddies  are  less  than 
satisfactory.  One  would  anticipate  great  difficulty  in  the  direct 
simulation  of  the  Navier-Stokes  equation,  since  in  the  past  150  years, 
only  about  75  analytic  solutions  have  been  found — a  tribute  to  the 
intractability  of  simultaneous  nonlinear  partial  differential  equations. 
The  pursuit  of  numerical  solutions  to  the  full  set  of  equations  is 
hindered  by  the  rapid  increases  in  degrees  of  freedom  (nodal  points  in 
finite  difference  methods)  with  Reynolds  number.  Schumann  et  al. 
[1980],   Liepmann   [1979],   and  Orszag   [1977]   have  all  noted  that  the 

9/4 
required  degrees  of  freedom  scale  with  Re    .   If  it  were  possible  to 

deal  with  a  flow  where  Re=100  with  1000  interior  mesh  points,  the  scal- 

4  7 

ing  law  suggests  that  Re=10  would  require  about  3  x  10  mesh  points. 

In  fact,  Orszag  notes  that  an  order  of  magnitude  increase  in  computa- 
tional power  will  permit  an  increase  in  Re  of  only  2.15  times  and 
Schumann  et  al.  note  that  a  brute-force  simulation  on  a  10  MIPS  machine 

4 
for  Re=10   would  require  about  3  years  computing  time.   There  are 

grounds  for  optimism,  however,  for  free  flows  and  flows  in  decay,  where 


either  the  Re  is  small  or  the  large  scale  motions  can  be  treated  inde- 
pendently of  the  dissipative  structure.  There  is  less  reason  to  be 
optimistic  about  success  in  direct  simulation  of  flows  about  objects; 
for  turbulent  boundary  layers  or  laminar  boundary  layers  undergoing 
transition,  the  rapid  changes  in  the  streamwise  direction  cause  dif- 
ficulties in  resolution.  Furthermore,  the  downstream  or  outflow 
boundary  conditions  would  appear  to  require  complete  specification  in 
the  wake  region — with  the  concommitant  danger  that  the  flow  will  be 
over-specified  and  the  results  of  the  direct  simulation  set  in  advance. 

Therefore,  there  is  need  for  numerical  methods  in  which  there  is 
not  a  rapid  increase  in  the  degrees  of  freedom  with  increasing  Reynolds 
number  for  the  full  set  of  Navier-Stokes  equations  and  which  also  does 
not  require  "numerical"  boundary  conditions  in  addition  to  those 
specified  by  the  original  problem.  Preferably,  it  should  be  possible  to 
apply  this  approach  easily  to  any  problem  of  interest. 


CHAPTER  II 
REVIEW  OF  LITERATURE 
For  the  above  reasons  the  decision  was  made  to  explore  the  spectral 
methods  described  in  the  monograph  of  Gottlieb  and  Orszag  [1977]  and  the 
review  of  Gottlieb  et  al  [1984a].  Spectral  methods  have  become  incre- 
asingly popular  in  recent  years  resulting  in  studies  of  transitional 
fluid  flows  (Kleiser  [1982],  Marcus  et  al.  [1982]  and  Orszag  and  Kells 
[1980]);  incompressible  fluid  flow  simulations  (Haidvogel  et  al.  [1980], 
Kumar  and  Yajnik  [1980],  Leonard  and  Wrey  [1982],  Noin  [1982],  Moin  and 
Kim  [1980],  Orszag  [1971,  1980]  and  Taylor  and  Murdock  [1980,1981]);  and 
most  recently  for  complicated  compressible  flow  fields  which  include 

shock  waves  (Gottlieb  et  al.  [1981,  1984b],  Hussaini  and  Zang  [1984], 
Sakell  [1984]  and  Salas  et  al.  [1982]).  Interest  in  spectral  methods  is 
due  in  part  to  the  increased  accuracy  available  for  a  given  number  of 
independent  degrees  of  freedom  (i.e.  the  number  of  mesh  nodes  for  a 
finite  difference  method)  in  comparison  to  finite  difference  techniques; 
for  example  see  Haidvogel  et  al.  [1980]  or  Dennis  and  Quartapelle  [1982] 
for  comparisons  of  accuracy.  The  spectral  methods  are  often  referred  to 
as  being  infinite-order  accurate,  that  is,  if  after  N  terms  the  error 
decreases  more  rapidly  than  any  power  of  1/N.  Orszag  [1971]  states  that 
it  is  possible  to  decrease  the  computational  time  and  storage  by  an 
order  of  magnitude.  The  advent  of  the  fast  Fourier  transform  in  the  mid 
60* s  facilitated  the  development  of  spectral  methods,  since  rapid 
evaluation  of  trigometric  polynomial  coefficients  became  possible.  Much 
of  the  work  done  with  respect  to  fluid  dynamics  has  involved  the  use  of 


the  pseudospectral  techniques  also  known  as  "collocation"  or  the  "method 
of  selected  points"  (Lanczos  [1956]). 

Spectral  methods  represent  the  solution  to  the  problem  as  a 
truncated  series  of  eigenf unctions  for  the  independent  variable.  The 
type  of  series  expansion  used  is  based  upon  the  type  of  boundary  condi- 
tion that  is  to  be  satisfied,  for  example,  a  periodic  boundary  condition 
would  suggest  the  use  of  a  Fourier  series  representation  for  the 
solution.  Gottlieb  and  Orszag  [1977]  and  Orszag  [1977]  suggest  the 
following  series  representations: 

Boundary  Conditions        Series  Representation 

periodic  exp(inx) 

inviscid  sin(nx)  or  cos(nx) 

no-slip  Chebyshev  (T  (x))  or  Legendre  (P  (x)) 

The  most  popular  series  representations  are  the  sine  series  and  the 
Chebyshev  series  due  in  part  of  the  ability  to  use  the  FFT  to  evaluate 
coefficients  or  their  inverse  transform.  Additional  series  repre- 
sentations have  been  used  for  special  problems;  for  example  Tang  [1979] 
used  surface  harmonics  as  the  expansion  function  for  two  dimensional 
flows  on  the  surface  of  a  sphere  and  Kasahara  [1977]  used  Hough  har- 
monics to  solve  shallow-water  equations  over  a  sphere.  Typically  the 
series  representation  was  an  eigenfunction  expansion  for  the  governing 
equation.  Gottlieb  and  Orszag  [1977]  discuss  the  convergence  properties 
of  various  series  expansions  and  conclude  that  for  "non-slip"  boundary 
conditions  that  Chebyshev  or  Legendre  series  provide  the  greatest  ac- 
curacy for  the  fewest  terms  in  the  expansion. 


Spectral  methods  can  be  classified  as  Galerkin,  Tau,  or 
Pseudospectral.  A  Galerkin  method  seeks  the  approximation  of  the  de- 
pendent variable,  u„(x,t),  in  the  form  of  a  truncated  series. 

N 

N 
uN(x,t)  =  Z  an(t)  <J>n(x)  (4) 

n=l 

where  4>  (x)  are  linearly  independent  functions  choosen  such  that  uN(x,t) 
satisfies  all  boundary  conditions.  The  coefficients,  a  (t),  are  deter- 
mined by  the  Galerkin  equations: 

5t      V  XVNdX>   "  d'  WX  +  Df  VdX  C»  -  1.    « ■) 

(5) 

where  L  is  a  linear  differential  operator.  The  tau  method  developed  by 

Lanczos  [1956]  employs  expansion  functions,  4>   ,   assumed  to  be  elements 

n 

of  a  complete  set  of  orthonormal  functions.  The  approximation  to  the 

solution,  u„(x,t)  is  expressed 

N 


uN(x,t)  =  Z  an(t)  $n(x)  (6) 


N+k 

Z 
n=l 

where  k  is  the  number  of  independent  boundary  constraints  that  must  be 


applied.   The  coefficients  a„  ,   to  a„  .   are  chosen  such  that  the  k 

N+l       N+k 

boundary  conditions  are  satisfied.  Whereas  the  expansion  coefficients 

a   to  a   are  the  approximate  solutions  for  the  differential  equation 

being  solved.   Collocation  methods  represent  the  solution,  u  ,  as  a 

N 

series  expansion: 


N 

u„  =  Z  a  <f>   (x)  (7) 

N     «   n^n 
n=l 

where  the  expansion  coefficients  a  are  the  solutions  of  the  equations 

Z  a  *  (x.)  =  u(x.)  (8) 

n  n  1       i 

Notice  that  the  coefficients  a  are  dependent  on  both  <t>     and  the  points 

n  n 

x   for  n=l ,   2 N.   Collocation  methods  can  be  applied  in  either 

n 

normal  or  spectral   space  (i.e.  solve  for  u  or  a  ,  respectively)  but, 

most  often  are  applied  in  normal  space  because  of  the  nonlinear  nature 
of  the  governing  equation.  The  greatest  difference  between  these 
methods  is  the  manner  in  which  boundary  conditions  are  handled  (see 
section  2  of  Gottlieb  and  Orszag  [1977]). 

Pseudospectral  (collocation)  techniques  are  the  most  likely  to  be 
appropriate  for  typical  governing  equations  of  transport  phenomena.  The 
principal  concept  behind  pseudospectral  calculations  as  stated  by  Orszag 
[1980]   is  to  simply  transform  freely  between  physical  (x.)  and  spectral 

(a  )  representations,   evaluating  each  term  in  whatever  representation 

that  term  is  most  accurately,  and  simply  evaluated.  Pseudospectral 
computations  have  several  advantages  over  spectral  algorithms:  1)  For 
complex  geometry  the  solution  of  the  spectral  (Galerkin)  method  requires 
at  least  twice  the  number  of  fast  Fourier  transforms  than  that  of  the 
pseudospectral  method  using  trigometric  series  representation;  2) 
Pseudospectral  techniques  have  significant  advantages  over  spectral 
techniques  when  solving  nonlinear  partial  differential  equations  espe- 
cially when  considering  computational  expense. 
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In  order  to  demonstrate  the  spectral  method  consider  the  one  dimen- 
sional transient  equation: 

9u   32u 

3t   .  2  (9) 

3x 

with  the  boundary  and  initial  conditions 

u(0,t)  =  u(/r,t)=0 

u(x,0)  =  f(x)  (10) 

the  analytic  solution  to  this  problem  can  be  shown  to  be 

00        Jt 

u(x,t)  =  Z  {  |  f  f (x)sin(nx)dx  }  sin(nx)exp(-n2t)      (11) 
n=l   n     0 

Therefore,   let  u  be  the  spectral  approximation  to  u,  so  that  the  solu- 
tion can  be  approximated  as 

N 
uN(x,t)  =  Z  a  (t)sin(nx)  (12) 

n=l 

Substituting  back  into  equation   (12)  and  using  the  Fourier  transform 
results  in 

dan(t) 

-dT~  =  _n  an(t)  <13> 

and  this  set  of  ordinary  differential  equations  is  solved  with  respect 

to  the  initial  conditions  a(0)  =  2/ir  J   f(x)sin(nx)dx  (n=l N) . 

0 

Gottlieb  and  Orszag  [1977]  have  shown  for  this  example  that  u(x,t)-u 

N 

(x,t)  goes  to  zero  more  rapidly  than  exp(-N  t)  for  any  t>0  and  N— >co. 
It  can  be  seen  that  the  spectral  approximation  of  u(x,t)  that  results  is 
a  truncation  of  the  exact  solution  to  N  terms.   In  contrast  to  this 


9 


exceedingly  good  approximation  it  i3  possible  to  construct  a  spectral 
method  that  will  arrive  at  extremely  poor  results;  therefore  care  must 
be  exercised  in  the  construction  of  a  spectral  approximation. 

The  pseudospectral  method  is  typically  chosen  over  other  spectral 
methods  when  considering  a  nonlinear  equation.   For  example  consider 

2 

3u      .   ,3u   9  u  ,   , 

aT  =  exP<xu)lx-  +  72  <14> 

3x 
let  the  approximation  u„(x,t)  to  u(x,t)  be 

N 

N 

uN(x,t)  =  Z  an(t)*n(x)  (15) 

n=l 

where  ♦  (x)  is  an  orthonormal  function.   Then  the  spectral  approximation 

will  be 

daQ(t) 

-St"  =  J  Ve^[atSVn1SVA  +  ^Vn)  dx  (16) 

It  is  seen  that  these  equations  for  a  (t)  are  computationally  complex 

due   to   the   resulting   integro-differential  equation  for  a  (t). 

n 

Therefore,   consider  the  pseudospectral  method  in  which,  N  collocation 

points   (Xj.Xg. . . .Xjj)   lying   within  the  computational  domain  are 

introduced.   The  approximation  (15)  is  forced  to  satisfy  the  governing 

partial  differential  equation  (14).   For  example,  the  following  steps 

would  be  followed;  1)   Determine  N  coefficients  a  (0)  such  that 

n 

N 
uN(x  0)  =  E  ■  (0)*  (x  )  (17) 

n=l         J 


10 


2)  Evaluate  each  term  of  the  governing  partial  differential  equation  in 
either  physical  or  spectral  space,  whichever  gives  the  most  accurate  and 
easily  obtainable  approximation.   For  example,   exp  (x  u  (x.,t))  is 

evaluated  in  physical  space  since  the  value  for  u„(x.,t)  is  known  and 

N   J 

the  partial  derivatives  are  evaluated  in  spectral  space  since  this 

results  in  the  most  accurate  representation.   3)   Integrate  in  time  with 

respect  to  uN(x.,t)   using  the  "leap-frog"  method  or  another  suitable 

choice.  4)  Repeat  steps  l)-3)  until  the  time  integration  is  completed. 
From  this  example  pseudospectral  methods  are  seen  to  be  much  easier  to 
apply  to  nonlinear  equations  than  corresponding  spectral  techniques. 

The  application  of  boundary  conditions  in  a  spectral  method  can 
determine  the  solution's  stability.  Gottlieb  et  al.  [1981]  state  that 
incorrect  boundary  treatments  may  give  strong  instabilities  in  contrast 
to  finite-difference  methods  in  which  this  would  appear  as  relatively 
weak  oscillations.  Moin  and  Kim  [1980]  note  that  fully  explicit  pseudo- 
spectral  solution  of  the  incompressible  Navier-Stokes  equation  have  an 
inherent  numerical  problem  for  viscous  flows  involving  solid  boundaries, 
due  to  enforcing  no-slip  conditions  at  the  walls.  Rudy  and  Strikwerda 
[1981]  presented  a  study  of  inflow  and  outflow  boundary  conditions  for 
compressible  Navier-Stokes  equations  of  flow  past  a  flat  plate.  This 
study,  conducted  for  finite-difference  methods  of  solution,  indicates 
that  errors  in  the  data  specified  at  the  inflow  boundary  condition  can 
significantly  affect  the  solution  obtained.  Gottlieb  and  Orszag  [1977] 
note  spectral  methods  are  extremely  sensitive  to  the  formulation  of 
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boundary  conditions;  for  example,  when  improper  boundary  conditions  are 
imposed,  the  solution  is  likely  to  be  "explosively"  unstable. 

Convergence  of  the  spectral  (Galerkin)  method  for  Navier-Stokes 
equations  has  been  demonstrated  for  a  Fourier  representation  with  peri- 
odic boundary  conditions  by  Hald  [1981].  Maday  and  Quarteroni  [1982] 
have  provided  stability  results  and  "optimal"  convergence  rates  for 
Galerkin  and  pseudospectral  approximations  (using  trigonometric  polyno- 
mials) for  the  stationary  Navier-Stokes  equation  with  periodic  boundary 
conditions.  Canuto  [1984]  analyzed  explicit  and  implicit  methods  of 
imposing  boundary  conditions  for  Chebyshev  and  Legendre  approximations 
of  elliptic  problems  ensuring  stability  and  convergence  for  these 
methods.  Canuto  and  Quarteroni  [1984]  give  stability  and  "optimal" 
convergence  rates  for  Chebyshev  collocation  approximations  of  variable 
coefficient  elliptic  problems  with  Dirichlet  or  Neumann-type  boundary 
conditions.  Gottlieb  [1981b]  has  demonstrated  stability  of  Chebyshev- 
pseudospectral  representations  for  parabolic  and  hyperbolic  equations 
with  variable  coefficients.  Gottlieb  and  Orszag  [1977]  discuss  al- 
gebraic stability  criteria  as  applied  to  spectral  methods.  Pasciak 
[1980]  investigated  spectral  and  pseudospectral  representations  of 
advection  equations  with  the  intention  of  introducing  a  framework  in 
which  finite  elements  analysis  can  be  applied  to  spectral  methods.  In 
addition  error  estimates  are  given  for  fully  discrete  explicit  pseudo- 
spectral as  well  as  semidiscrete  spectral  and  pseudospectral  methods. 

Some  additional  considerations  warrant  attention.  In  the  case  of  a 
discontinuity,  the  rate  of  convergence  in  the  region  of  the  discon- 
tinuity is  seriously  degraded,   but  spectral  approximations  are  still 
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normally  more  accurate  than  corresponding  finite  difference  repre- 
sentations (Gottlieb  and  Orszag  [1977]).  Additionally,  the  error  is 
localized  better  by  the  spectral  method  such  that  less  local  dissipation 
is  required  to  smooth  the  discontinuities.  If  dissipation  or  a  filter- 
ing technique  is  not  used  then  Gibb's  phenomenon  is  seen  in  the  region 
of  the  discontinuity,  and,  the  resulting  error  from  the  lack  of  smooth- 
ing pollutes  the  solution  globally  (Osher  [1984]).  Majda  et  al.  [1978] 
have  shown  for  general  linear  hyperbolic  Cauchy  problems  with  nonsmooth 
initial  data  that  the  appropriate  smoothing  techniques  applied  to  the 
equation  results  in  stability  and  that  this  smoothing  combined  with 
smoothing  of  the  initial  data  gives  rise  to  infinite  order  accuracy  away 
from  the  discontinuities  of  the  exact  solution.  The  most  prominent 
examples  of  discontinuities  that  have  been  examined  by  spectral  methods 
are  shock  waves  (Gottlieb  et  al.  [1981,1984b],  Hussaini  and  Zang  [1984], 
Sakell  [1984],  Salas  [1982],  Streett  et  al.  [1985],  and  Taylor  et  al. 
[1981]).  A  variety  of  filters  and  dissipation  functions  have  been  used 
to  damp  the  oscillations  occurring  in  the  solution  due  to  the 
discontinuity.  For  example,  one  sided  Schumann  filter  (Gottlieb  et  al. 
[1981]),  von  Hann  window  filter  (Hussaini  and  Zang  [1984]  and  Salas 
[1982]),  second  and  fourth-order  artificial  viscosity  (Sakell  [1984] 
and  Streett  et  al  [1985]),  artificial  density  (Street  et  al.  [1985])  as 
well  as  the  method  of  Boris  and  Book  [1976],  orginally  developed  for  the 
construction  of  finite  difference  algorithms  involving  strong  shocks 
(Taylor  et  al.  [1981]),  have  produced  accurate,  smooth  solutions  for 
this  type  of  problem.  In  addition,  Gottlieb  et  al.  [1981]  used  a  low- 
pass  spectral  filter  to  remove  the  high  frequency  waves  that  lead  to 
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instability  and  a  "cosmetic"   filter  (a  one-sided  Schumann  filter)  to 

find  a  nonoscillatory  numerical  solution.   Gottlieb  [1985]  provides  a 

brief  review  of  recent  advancements  in  the  field  of  compressible  flow 

problems.   Streett  and  Bradley  [1985]  discuss  briefly  some  applications 

of  spectral  methods  in  aerodynamics. 

Symmetric  flow  past  a  flat  plate  has  been  treated  by  Orszag  [1971] 

and  Taylor  and  Murdock  [1980] , [1981] .   Planar  flow  has  been  investigated 

by  Moin  and  Kim   [1980],   Orszag  and  Kells  [1980]  and  Kleiser  [1982]. 

Taylor  and  Murdock  studied  the  flow  over  a  flat  plate  in  a  range  of 

5         5 
1.2x10  <Re<3.8xlO  ,  which  is  well  within  the  turbulent  flow  regime.   The 

velocity  profile  was  solved  as  perturbation  about  the  Blasius  velocity 
profile.  Taylor  and  Murdock  used  a  mesh  of  17x17  to  solve  for  the 
velocity  profile,  while  Orszag  used  nine  collocation  points  to  ap- 
proximate a  one  dimensional  boundary  layer,  both  studies  resulted  in 
accurate  approximations. 

The  treatments  of  transonic  flow  by  Gottlieb  et  al.  [1984b]  and 
Streett  et  al.  [1985]  are  interesting  in  that  spectral  methods 
(Chebyshev  collocation)  are  being  used.  The  treatment  of  transonic  flow 
over  an  airfoil  can  be  viewed  as  state-of-the-art  application  of 
spectral  methods  because  of  the  difficulty  posed  by  the  sharp  shock 
gradients  and  the  computational  competition  with  finite  difference 
methods.  Again,  it  is  seen  that  separation  causes  difficulty  as  Streett 
et  al.  note  that  the  potential  equation  which  is  being  solved  becomes  a 
mixed  elliptic-hyperbolic   type   and  admits  weak   solutions   with 
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discontinuities.  Streett  et  al.  use  artificial  viscosity  with  a  direc- 
tional bias  introduced  in  the  potential  equation  in  the  supersonic 
region  to  suppress  the  appearance  of  compression  and  expansion  shocks 
due  to  the  presence  of  a  supersonic  bubble.  It  is  apparent  that  the  use 
of  nultigrid  techniques  has  made  spectral  methods  for  steady  compres- 
sible flow  competitive  with  finite  difference  methods  for  problems  of 
aerodynamic  interest. 

A  very  thorough  review  of  application  and  major  developments  of 
spectral  methods  from  approximately  1977  to  1983  is  provided  by  Gottlieb 
et  al.  [1984a].  This  review  presents  all  major  aspects  of  applying 
spectral  methods  to  Navier-Stokes  equations.  Taylor  et  al.  [1985]  have 
proposed  a  method  for  solving  the  primative  three  dimensional  (3-D) 
Navier-Stokes  equations  without  introducing  the  Poisson  equation  for 
pressure.  This  was  accomplished  by  integrating  the  equation  of  con- 
tinuity for  the  normal  velocity  component  with  respect  to  the  surface 
and  integrating  the  corresponding  equation  of  motion  to  evaluate  the 
pressure  term.  This  procedure,  applied  to  evaluate  boundary  layer 
stability,  was  found  to  be  the  same  as  a  marching  solution  of  a  Poisson 
equation  and  therefore  unstable.  A  second  method  was  proposed  in  which 
the  normal  velocity  component  was  evaluated  from  both  the  equation  of 

continuity  (V  )  and  the  equation  of  motion  (V  )  and  an  iterative  method 
c  n 

of  the  following  form: 

p-i ,  p» .  a  -a,vv  (18) 

was  used  to  evaluate  the  pressure  at  the  new  time.   They  found  this 
approach  to  converge  very  slowly  due  to  its  explicit  nature.   Malik  et 
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al.  [1985]  describe  a  Fourier-Chebyshev  spectral  method  for  the  incom- 
pressible Navier-Stokes  equations.  The  algorithm,  combines  a  fully 
spectral  scheme  with  a  Fourier-finite  difference  method  in  evaluating 
most  wall  bounded  shear  flows.  Zang  and  Hussaini  [1985]  have  employed 
this  algorithm,  with  appropriate  modifications,  using  Fourier-Legendre 
expansions  to  model  incompressible  channel  flow.  Satisfactory  results 
were  obtained.  It  should  be  noted  that  while  Legendre  series  and 
Chebyshev  series  require  approximately  the  same  number  of  terms  to 
converge  to  the  solution,  (Gottlieb  and  Orszag  [1977])  Legendre  series 
are  rarely  used  in  practice.  While  the  FFT  is  much  faster  than  matrix 
inversion  for  large  matrices  as  the  matrix  size  decreases  the  efficiency 
of  the  FFT  also  decreases  until  it  is  approximately  that  of  matrix 
inversion;  therefore  Legendre  representations  should  be  competitive 
Chebyshev  representations  when  "few"  terms  are  required  for  an  accurate 
approximation.  Taylor  [1984]  has  shown  that  for  N  <  64  that  Crout's 
method  is  equivalent  or  better  than  the  FFT.  Also  Zang  and  Hussaini 
note  that  evaluation  of  derivatives  by  matrix  multiplication  (Legendre 
collocation  matrices)  is  faster  for  N=32  than  machine  language  Chebyshev 
transforms . 

Many  methods  are  being  used  to  solve  spectral  representations,  for 
example,  Deville  et  al.  [1982]  examined  time  integration  for  the  non- 
linear Burger  equation.  They  conclude  that  only  two  methods  are 
practical  choices:  The  first  employs  a  second  order  Adams-Bashforth 
method  to  integrate  the  non-linear  terms  and  a  second  order  Crank- 
Nicolson  scheme  to  integrate  the  linear  terms  and  the  second  method 
utilizes  predictor  and  corrector  steps  and  both  finite  difference  and 
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pseudospectral  methods.  Taylor  [1984]  outlines  a  finite  difference 
predictor  and  spectral  corrector  approach  for  integration.  An  ADI 
method  has  been  used  to  solve  a  Poisson  equation  (Haidvogel  and  Zang 
[1979]).  Haidvogel  and  Zang  described  a  second  algorithm  for  Poisson's 
equation  in  which  the  tau  method  is  used  and  the  resulting  coefficient 
matrix  for  the  second  order  derivative  is  diagnolized.  This  approach  is 
better  suited  to  time-dependent  problems  and  is  an  order  of  magnitude 
more  efficient  than  the  ADI  algorithm.  Haldenwang  et  al.  [1984]  have 
used  this  algorithm  for  a  3-D  Helmholtz  equation.  The  use  of  non- 
homogeneous  boundary  conditions  with  the  resulting  equations  for 
evaluation  of  derivatives  is  discussed  in  Appendix  I  of  the  Haldenwang 
et  al  paper.  Sharp  and  Harris  [1984]  have  investigated  the  combination 
of  pseudospectral  methods  with  parametric  differentiation  as  applied  to 
the  Falkner-Sakan  equation.  From  CPU  times  given  for  evaluation  of  the 
equation,  it  does  not  appear  that  the  method  of  parametric  differentia- 
tion combined  with  a  pseudospectral  expansion  is  competitive  with  the 
pseudospectral  method.  Zebib  [1983]  has  proposed  a  new  Chebyshev  method 
which  represents  the  highest  derivative  as  a  Chebyshev  series  expansion 
with  lower  derivatives  being  obtained  by  integration.  This  technique 
can  be  viewed  as  being  related  to  the  tau  method  in  that  additional 
terms  arising  from  the  integration  are  used  to  satisfy  boundary 
conditions.  From  the  examples  presented  the  method  appears  to  be  fully 
as  accurate  as  more  traditional  Chebyshev  techniques. 

For  cases  where  fully  spectral  techniques  are  difficult  to  apply  or 
inappropriate  for  the  problem  domain  other  avenues  are  available  for 
their  application.   One  possible  method  is  the  combination  of  finite 
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difference  and  pseudospectral  methods  for  separate  directions.  Reddy 
[1983]  used  such  an  approach  in  solving  3-D  flowfields  over  missle 
shaped  configurations  at  moderate  angles  of  attack.  In  this  approach  a 
Fourier  sine  series  was  used  in  the  circumferential  direction  with 
finite  difference  representation  for  the  other  directions.  Reddy  es- 
timated that  a  40%  increase  in  efficiency  could  be  obtained  for  a 
relatively  sparse  grid  with  respect  to  the  spectral  representation. 
Metivet  and  Morchoisne  [1982]  have  investigated  splitting  the  domain  on 
which  the  problem  is  defined  into  a  finite  number  of  problems  on 
subdomains.  Metivet  and  Morchoisne  note  the  main  considerations  of  this 
approach  are  "patching"  (i.e.  ensuring  that  derivatives  at  the  bound- 
aries of  adjacent  subdomains  are  continuous)  and  modifications  of 
equations  due  to  mapping  of  the  subdomains.  The  proposed  method  has  the 
flexibility  of  finite  element  and  finite  difference  methods,  but  main- 
tains the  accuracy  of  monodomain  spectral  techniques. 

It  is  evident  that  spectral  methods  can  be  applied  competitively  to 
many  aerodynamic  and  hydrodynamic  flow  simulations.  Since  equations  for 
heat  and  mass  transfer  are  often  similar  to  those  of  fluid  flow,  it  is 
likely  that  spectral  methods  can  be  applied  profitably  in  those  cases  as 
well.  The  purpose  of  this  work  is  to  investigate  the  application  of 
spectral  methods  to  problems  in  transport  phenomena  and  to  assess  the 
suitability  and  computational  efficiency  of  such  methods  when  applied  to 
important  nonlinear  equations. 
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CHAPTER  III.   APPLICATION  OF  THE  CHEBYSHEV  PSEUDOSPECTRAL  METHOD 
In  order  to  discuss  how  Chebyshev  collocation  is  applied,  it  is 
advantageous  to  have  a  framework  for  the  evaluation  of  coefficients  and 
derivatives.   A  Chebyshev  series  is  defined  on  the  interval  (-1,1)  with 
the  interpolation  points  most  commonly  chosen  as 

Xj  =  cos  ^j      (j  =  0,1 N)  (19) 

for  a  Nth-order  representation.  A  Chebyshev  polynomial  of  degree  n  is 
defined 

T  (x)  =  cos(n  cos   (x) ) 

therefore  it  follows 

Tn(Xj)  =  cos  ^      (j  =  n  =  0,1 N)  (20) 

Based  on  these  definitions  the  approximation  to  the  solution  U„(x.) 

Nv  J' 

is: 


n=0 
where   a   are  the  polynomial  coefficients.   The  coefficients  are 

evaluated  by  application  of  a  Fourier  transform  resulting  in: 

N 
Cnan  =  N  Z  Cj  W  cos  ^       (n=0,1 N)      (22) 

where  cQ  =  cN  =  2  and  c  =  1  for  0<j<N.   This  transform  can  be  expressed 
in  terms  of  matrices 
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a  =  T     U„ 


(23) 


where 


a  = 


"1 


UN  = 


UN   <V 

uN  (xl} 


UN    (XN 
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T     =  — 
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Now  consider  the  derivative  of  U  which  can  be  expressed  in  terms  of 
polynomial  series  the  following  manner 

5x"<Vx»  =  Vn5x-<VX» 
n=0 


for  a  Chebyshev  representation  this  can  be  rewritten 

h   <VX»  "  \   bnVx) 
n=0 


(24) 


where  b   =  0.  The  relation  between  the  coefficients  a  and  b  can  be 
w  .  n     n 

shown  to  be 
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2  N 

b  =  —  Z 

n   c  , 

n  p=n+l 


P  a. 


(25) 


p+n  odd 
based  on  the  recurrence  property 


2T  (x)  = 
n 


c  T   . '(x)   w  T    "(x) 
n  n+1  *  '    n  n-1  v  ' 


n+1 


n-1 


(26) 


where  c  =w  =0  if  n<0,   c_=2,  w  =1,  and  c  =w  =1  for  n>0,  and  the  prime 
n  n  u     o        n  n 

indicates  differentiation.   The  coefficient  representation  for  the 

second  derivative  can  be  shown  to  be: 


dn  =  c" 


N 
Z 
n   p=n+2 
p+n  even 


,  2   2, 
p(p  -n  )at 


(27) 


where 


2         N 

7T  <UN>  -  *     dn  Tn^) 
dx        n=0 


(28) 


and  dN=dN-l=  °"   The  evaluatlon  of  derivative  coefficients  can  be 

represented  in  matrix  fashion.   Consider  the  coefficients  of  the  first 
derivative 


where 


b  -  D<»  . 


r  b 


(29) 


and 


N  J 


dJ.'  =0  if  i>j  or  i+j  even 
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otherwise 


DJ}*  =  ?!  for  i=j=0,l N 

i  j    c . 

The  second  derivative  coefficients  can  be  expressed 

d=D(1)b=  tD(1)D(1)a]  =D(2)a 
where 


(30) 


d  = 


r  d 


LdNJ 


It  can  be  seen  that  any  order  derivative  can  be  evaluated  by  raising 

D     the  appropriate  power,   q,   where  q  is  less  than  N.   Ku  and 
Hatziavranidis  [1985]  have  taken  this  approach  a  step  further  by  noting 


_n_(l)_s  . 
=  T  D   T  U 


N 


(31) 


or  for  a  qth  order  derivative 

UM(q)  -T"  (D(1))qT8  UM 


(32) 


where 
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VV      VV      T2(xo>    •  •  •  w 


W      Tl(Xl)      T2(Xl)    .  .  .   TN(Xl) 


t0(X2)        yX2)        t2(x2)     . 


TN(X2) 


[T^)      Ti(Xn)      T2(Xn)    .  .  .   Tn(Xn) 


They  note  that  this  procedure  is  as  accurate  as  the  FFT  in  evaluating 
derivatives  and  is  more  efficient  in  terms  of  computational  time  and 
implementation.  The  basic  framework  for  evaluation  of  derivatives  and 
coefficients  has  been  established. 

There  are  limitations  to  the  above  approach.  One  of  the  most 
notable  exceptions  is  the  application  of  Robin's  boundary  conditions  at 
the  interface  of  two  calculated  domains  (i.e.  flux  equated  across  the 
boundary) .  Often  the  derivatives  of  the  dependent  variable  for  the 
first  one  or  two  points  beyond  the  boundary  are  not  very  accurate  in  the 
pseudospectral  representation  due  to  aliasing.  For  a  cosine  function 
cos(ux),  the  frequencies  u  and  -<j  are  indistinguishable  and  are  said  to 
be  aliases  of  each  other.  More  specifically  when  a  Fourier  transform  is 
used  to  transform  a  function  into  spectral  space  and  back  again  the 
result  is: 

w(k)  =  w(k)  +  w(k+2K)  +  w(k-2K)    (IlklKK)  (33) 

where  w(k)   is  the  aliased  representation  of  w(k)  and  N  =  2K.   The  last 
two  terms  arise  from  the  fact  that  exp[i(k±N)x  ]  =  exp[ikx.]  for  all 

J  J 
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integral  j  and  k.  It  should  be  noted  that  for  IlklKK  only  one  of  the 
two  terms  due  to  aliasing  can  be  non  zero.  Therefore,  the  values  of  the 
derivative  at  the  boundary  must  be  forced  to  satisfy  the  required 
condition.  This  can  be  accomplished  either  by  direct  substitution  of  the 
derivative  value,  if  possible,  or  by  the  use  of  the  tau  method. 
Additional  summation  terms  are  needed  for  the  cases  where  aliasing  is 
significant  and  hence  has  to  be  removed.  A  matrix  approach  may  be 
difficult  to  apply  for  the  case  of  nonlinear  problems.  Therefore,  in 
certain  cases,  it  may  be  easier  to  eliminate  the  aliasing  if  the 
governing  equation  is  represented  in  terms  of  the  series  coefficents. 
Another  example  is  when  a  spectral  filter  is  desired  to  eliminate 
"noise"  from  the  solution  arising  from  a  discontinuity.  Again,  it  is 
easier  to  use  the  series  coefficents  in  a  spectral  filtering  method  than 
to  use  the  dependent  variable  values. 

Using  the  preceeding  approach  classical  problems  in  fluid  and 
thermal  transport  and  two  additional  problems  dealing  with  diffusion 
were  studied.  The  classical  problems  considered  were;  1)  Start-up 
laminar  velocity  profile  in  a  circular  tube;  2)  Transient  temperature 
profile  of  a  wire  filament  electrically  heated;  3)  The  classical 
Graetz  problem.  Diffusion  in  which  diffusivity  was  a  function  of 
concentration  was  modeled  for  a  planar  geometry.  The  last  problem 
considered  was  transient  diffusion  and  reaction  within  a  spherical 
pellet,  with  nonlinear  Michael is-Menton  kinetics  (typically  seen  for 
many  enzymes ) . 
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Start-up  Laminar  Velocity  Profile 
The   simplest  problem  considered  is  the  start-up  of  the  laminar 
velocity  profile  in  a  circular  tube.   This  problem  is  considered  simple 
because  of  the  linearity  of  the  governing  equation,  which  is  given  by: 

3v   Po  ~  PL     1  8  .   8v,  ,   . 

p  -pr 

v  is  the  fluid  velocity,  p   is  the  density  of  the  fluid,  t  is  time, 

is  the  applied  pressure  gradient,  fi   is  the  fluid  viscosity  and  r  is  the 

tube  radius.   The  following  initial  (I.C.)  and  boundary  conditions 

(B.C. ' s)  were  used: 

I.C.       v(r.O)  =0;    0  <  r  <  R  (35) 

B.C. 's     v(R,t)  =0;    t  >  0  (36) 

8v 

—  (0,t)  =  0;  t  >  0  (37) 

The  analytic  solution  to  this  problem  expressed  in  terms  of  the 
dimensionless  varibles  <j>,   £  and  T  where 

a,  = y. .     * «  £.    T  .  tSL. 

T  o  '  S  ft  >      *         o 

(Po-pL)R':/4/iL         K        pR* 
is 

2       •  Jo  (an?)  6XP  [_Gtn  T] 

♦<?.t)  =  (lY)   -  8  Z  -2 S. 2 —  (38) 

n=l      o  J„(cr  ) 
n  lv  n' 

The  an  are  chosen  such  that  the  resulting  zeroth  order  Bessel  function 

(J0(an))   ls  zero.   The  pseudospectral  approximation,  <j>         to  the 

solution,  <f>,     was  expressed  by  a  Chebyshev  polynomial  expansion.   The 
governing  equation  was   integrated  with  an  Euler  predictor.   Table  1 
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illustrates  the  results  obtained  for  this  problem  with  nine  and 
seventeen  expansion  terms  and  t  =  0.05.  Table  2  illustrates  the  results 
obtained  for  this  problem  with  seventeen  and  twenty-five  expansion  terms 
and  T  =  0.10,  results  for  a  finite  difference  representation  are 
included. 

Tables  1  and  2  illustrate  that  fairly  accurate  approximations  are 
obtained  for  the  central  region  of  the  tube.  Table  2  indicates  that  no 
significant  increase  in  accuracy  is  obtained  by  increasing  the  number  of 
coefficents.  An  anomaly  seen  for  the  psuedospectral  approximations  is 
the  relative  inaccuracy  of  the  solution  near  the  boundary  in  relation  to 
that  of  the  central  region.  The  cause  of  this  inaccuracy  may  be  due  to 
the  integration  method  chosen  or  presence  of  a  discontinuity  at  the 
boundary  when  the  pressure  gradient  is  intitiated.  The  results 
presented  in  Table  2  suggest  the  integration  method  used  was  relatively 
accurate  and  therefore  the  problem  in  the  boundary  region  is  due  to  the 
initial  discontinuity.  Possible  methods  to  deal  with  the  problem  of  the 
initial  discontinuity  include  using  a  digital  filter  to  remove  any 
"noise"  that  may  occur  near  the  boundary,  the  use  of  a  more  accurate 
integration  scheme,  or  use  of  a  highly  accurate  finite  difference 
representation  to  initialize  the  values  of  the  pseudospectral 
approximation  at  some  point  beyond  the  initial  starting  values. 
Aliasing  is  not  a  factor  in  this  case,  since  as  Orszag  [1972]  has  noted, 
aliasing  terms  lead  to  the  numerical  solution  being  susceptible  to 
numerical  instability. 
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Table  1. 


Comparison  of  pseudospectral  predicted  unsteady  laminar  flow 
to  the  analytical  solution  with  four  terms  of  the  series 
being  retained,  T  =  0.05. 


1.0 


0.9239 


0.7071 


0.3827 


0.0 


4>    (0.05,?)     0.0 


0.1054 


0.2056 


0.1941 


0.2000 


4>8(0.05,£)     0.0 


*  Error 


0.07848 


0.0     25.5 


0.1796 


112.6 


0.1998 

0.2032 

3.0 

1.6 

0.1959 

0.1976 

0.93 

1.2 

<f>16(0.05,?)    0.0 


*  Error 


0.06040      0.1579 
0.0     42.7         23.2 
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Table  2.  Comparison  of  pseudospectral  and  finite  difference  predicted 
flow  to  the  analytical  solution  with  four  terms  of  the 
series  being  retained,  T  =  0.10. 


1.0 


0.9239 


0.7071 


0.3827 


0.0 


<*>(o.io,e) 


0.0 


0.1182 


0.2861 


0.3577     0.3851 


<fr16(0.10,£)    0.0 


0.08447 


0.2504 


0.3584 


0.3885 


%   Error 


0.0 


28.5 


12.5 


0.20 


0.88 


<fr24(0.10,£)    0.0 


0.08443 


0.2503 


0.3582 


0.3886 


%   Error 


0.0 


28.6 


12.5 


0.14 


0.91 


Finite 
Difference 
(0.01  radial 
spacing 


0.0 


0.1115 


0.2990 


0.3864 


0.3991 


%   Error 


0.0 


5.7 


4.5 


8.0 


3.6 
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Transient  Temperature  Profile  of  a  Wire  Filament 
Electrically  Heated 
The  problem  considered  is  the  transient  temperature  profile  of  a 
cylinderical  wire  filament  electrically  heated  with  only  radiative  heat 
losses  considered.    It  is  assumed  that  the  physical  properties  of  the 
wire  are  constant.  The  governing  equation  is 

|J  =  a  2_I  +  0  (T*  _  T4)  (39) 

3t     3x2       m 


where 


a   =  —  ;  thermal  diffusivity 


and 


^   c<Jo 

T   -  (l!fi   ♦  T4)* 
m   xupea    o 

k  is  the  thermal  conductivity,  c  is  the  specific  heat  capacity  of  the 

filament  material,  6     is  the  density,  p  is  the  periphery  of  the  cross 

section  of  the  filament,  6  is  the  total  emissivity  of  the  surface,  cr  is 

the  Stefan-Boltzmann  constant,  u     is  the  cross  section  area,  I  is  the 

heating  current,  p  is  the  specific  resistance  and  T  is  the  temperature 

of  the  chamber  walls  in  which  the  heating  is  occuring.   The  following 
I.C.  and  B.C.'s  were  assumed: 
I.C. 

T(x,0)  =  T    ;     -1  <  x  <  1  (40) 
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B.C.  's 

T(l,t)  =  T(-l,t)  =  Tq  ;   t>0  (41) 

An  analytic  solution  for  the  central  region  of  the  wire  is  given  by  Jain 

and  Krishnan  [1955] 

T  -T 

Y^  exp[-2   tan  (— )]    =  C  exp    [-  |-]                                                          (42) 

m  ma 


where 


and 


t   =  (wV1 

0  IB 


T  -T  T 

C  =  J-_2  eXp[-2  tan"1    (f)] 
mo  m 


The     system  was  modeled  for  a  wire  two  centimeters (cm)    in     length  with  a 

=     0.1,      0     =     8.6X10"11,      T        =      1182. 5°K,      T       =   281°K,    t      =   1.762  and 

m  oo 

C=0.3863.   The  temperature  profile  obtained  for  a  time  of  0.04  seconds 

is   displayed  in  Figure  1.   This  profile  was  obtained  for  both 

pseudospectral  and  finite  difference  representations  using  a  second 

order  Adams-Bashforth  predictor  with  time  increments  of  0.0001  seconds. 

The  values  obtained  from  the  numerical  solutions  for  the  central  region 

of  the  wire  agree  exactly  with  that  obtained  from  the  analytic  solution. 

The  finite  difference  representation,   with  nodes  spaced  0.01  cm  and 

using  symmetry  about  x=0,   required  approximately  1.17  seconds  per 

iteration   on  a  Zenith  Z-122  microcomputer  using  compiled  Basic. 

Approximately  0.42  seconds  per  iteration  was  required  for  a  twenty-five 

term  pseudospectral  representation  on  the  same  machine.   In  both  cases 
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Figure  1.  Tenperature  Profile  within  an  Electrically  Heated  Wire 
Filament  Obtained  with  Finite  Difference  and  Pseudospectral 
Methods  for  t=0.04  seconds. 
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the  number  of  terras  retained  for  the  physical  representation  may  be 
reduced  without  affecting  the  accuracy  of  the  numerical  solution,  thus 
increasing  the  speed  of  the  computation. 

The  Classical  Graetz  Problem 

The   classical  Graetz  problem  refers  to  the  analysis  of  laminar 

flow  heat  transfer  in  a  tube.   For  the  case  of  a  circular  tube  the 

following  conditions  are  assumed;    1)   The  heat  capacity,   thermal 

conductivity,  viscosity,  and  density  of  the  fluid  are  constant;   2)  The 

laminar  velocity  profile  is  established  before  entering  the  region  in 

which  heat  transfer  will  occur;   3)  At  x=0  the  wall  temperature  will 

change  from  that  of  the  uniform   fluid  temperature,  T  ,  to  the  new 

o 

value,   T  ,   and  is  constant  at  this  value  for  x>0.   The  governing 

equation  for  the  heat  transfer  occuring  in  a  cylindrical  tube  neglecting 
thermal  conduction  along  the  axis  is 
3T    k   ,1  3   .   3T., 

u  ai  =  c"p  [r  §7  (r  aF)]  <43> 

p 

where  u  is  the  velocity  at  the  radial  position  r,  k  is  the  thermal 
conductivity,   C  is  the  heat  capacity,  p   is  the  density  of  the  fluid,  T 

is  the  local  temperature  and  x  is  the  axial  postion.  The  laminar 
velocity  profile  is  defined  to  be 

u  =  2U  [1  -  (J-)2]  (44) 

w 

U  is  the  maximum  velocity,   rw  is  the  radial  distance  to  the  wall. 
Boundary  conditions  employed  are 
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T(r,0)  =  TQ;     0  <  r  <  rw  (45) 

T(rw,x)  =  Tw;    x  >  0  (46) 

3T(0,x)    „  ,      v 

— -£-?■   =  0     x  >  0  (47) 

The  solution  takes  the  form: 

2  (— ) 
T  -  T     09  -/8  vr  ' 

tH-T  '     \   Cn*n  <]H  -Pt-V1-^  <48> 

w    o   n=0       w 

the  values  for  c  ,  o>  and  /S  are  given  by  Knudsen  and  Katz  [1958]  for 

n   n      n 

n=0,l  and  2  which  is  sufficient  for  most  situations.   The  governing 

equation   was   integrated  using  a  second  order  Adams-Bashforth  method 

-5 
with  axial  increments  of  1x10   .   It  was  found  by  trial  and  error  that 

within  the  first  30  steps  larger  axial  increments  (e.g.  doubling  the 

increment  size)  resulted  in  numerical  instability.  A  digital  filter  of 

the  form 

f'(x)  =  (f(x+2h)  +  2f(x+h)  +  9f(x)  +  2  f(x-h)  +  f(x-2h))/15    (49) 

was  used  initially  to  damp  oscillations  occurring  from  the  presence  of  a 

discontinuity  at  the  entrance.  This  filter  was  modified  to  account  for 

the  varible  node  spacing  that  results  for  Chebyshev  representations 

f'(xi)  =  (VUi-2)  +  'l  f(Vl)  +  9f(xi)   +  2f(Xi+l) 

+  f(x.+2))/(12+*1+J2)  (50) 

where 

*i  =  2  Ui  "  Vi)/(xi+i  -  V 

and 
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S2   =  (Xi-2  '  Xi)/(Xi  "  W 
for  i>l.  For  i=l  the  following  filter  was  used 

f'Uj)  =  (0.25  f(XQ)  +  10  f(x1)  +  f(x  ))/11.25  (51) 

Weighting  was  chosen  such  that  the  boundary  condition  at  the  wall  was 
not  over  emphasized  by  the  filter.  The  filters  were  used  only  for  the 
first  seven  axial  increments.   Figure  2  presents  numerical  results  in 

relation  to  values  obtained  from  the  series  summation  with  x/r  =  0.075. 

w 

It  can  be  seen  that   good  agreement  is  obtained  considering  the 

existence  of  a  discontinuity  at  the  entrance  region.   More  accurate 

results  may  possibly  be  obtained  by  the  use  of  different  filtering 

methods,   different  axial  spacing  between  filterings  or  the  use  of  more 

accurate  integration  methods  once  away  from  the  entrance  region.   Ku  and 

Hatziavramidis   [1984]   examined  the  classical  Graetz  problem  using 

Chebyshev  psuedospectral  representation  in  the  radial  direction  and 

either  a  Chebyshev  finite  difference  or  Chebyshev  finite  element  method 

in  the  axial  direction.   The  use  of  Chebyshev  techniques  in  the  axial 

direction  was  accomplished  by  transforming  the  problem  from  an  infinite 

domain  on  the  x-axis  to  a  finite  domain.  The  transformation 

z  =  J  tan-1  (— )  (52) 

w 

mapped  the  problem  onto  a  domain  of  z  ■  t   0.5,  where  K   is  a  constant 

that  may  be  varied  at  will.   Both  the  Chebyshev  finite  difference  and 

finite  element  techniques  gave  excellent  results  in  general. 
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Figure  2.  Solutions  to  the  Classical  Graetz  Problem  from  the  Analytic 

Series  Summation  and  the  Pseudospectral  Method  for  x/r,  =  0.075. 
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Diffusion  in  a  Planar  Slab 
This  problem  considers  the  transient  diffusion  of  ethanol  in  a 
porous  planar  slab  with  concentration-dependent  diffusivity.   Values  for 
the  diffusivity  of  ethanol  in  water  (Bird  et.  al.  [I960])  were  fit  to  a 
third  order  power  series  of  the  form 

D  =  a  +  bXE  +  cXg  +  dXg  (53) 

where  a  =   1 .443935xl0~5 

b  =  -7.373196xl0~5 

c  =   1.615881xl0~4 

d  =  -7.973184xl0~5 

as   determined  with  a  statistical  software  package  (3-D  Designer 

Statistics  for  Micros)  where  D  is  diffusivity  and  X_  mole  fraction  of 

E 

ethanol.   The  governing  equation  for  transient  diffusion  of  ethanol  in 

a  slab  infinite  along  the  x  and  z  axes  is 

8XE   dD  9XE   na\  ,   t 

ir  =  d^"^     "17  (54) 

An  expression  relating  the  derivative  of  the  diffusivity  to  ethanol  mole 
fraction  can  be  introduced  through  the  chain  rule 

42  =  42-^  (55) 

dy   3XE  9y  (55) 

Therefore  the  governing  equation  can  be  expressed  in  terms  of  ethanol 
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mole  fraction  only 

3X  3X  9  X 

-g|  =  (b  +  2cXE  +  3dXg)  (g-5)2  +  (a  +  bX£  +  cX2  +  dx!)  §   (56) 

3y 

Initial  and  boundary  conditions   impressed  upon  a  slab  2  centimeters 

thick  are 

X£  (y,0)  =  1;  -Ky<l  (57) 

XE(+l,t)  =  XE(-l,t)  =  0;    t>0  (58) 

Both  a  second-order  Adams-Bashforth  and  a  fourth-order  Runge-Kutta 
scheme  were  used  to  integrate  the  governing  equation.  Chebyshev 
representations  of  17,  25  and  33  terms  were  used  to  evaluate  the 
solution  in  addition  to  a  finite  difference  method  used  for  comparison. 
A  second-order  Adams-Bashforth  method  was  initially  applied  with  digital 
filtering  to  evaluate  the  governing  equation.  This  particular  technique 
displayed  severe  discrepancies  from  the  finite  difference  solution  and 
even  among  the  various  pseudospectral  solutions.  Integrating  to  t=1000 
seconds  invaribly  led  to  the  first  term,  away  from  the  wall,  being 
approximately  0.62  ±  0.02  for  each  approximation.  Subsequently,  a 
fourth-order  Runge-Kutta  method  was  applied  to  assess  the  validity  of 
the  Adams-Bashforth  method  that  was  used.  This  approach  gave  the  same 
results  as  the  Adams-Bashforth  method  to  within  5*.  Integrating  to 
t=5000  seconds,  with  17  terms,  still  resulted  in  the  term  next  to  the 
wall  being  approximately  0.62  ±  0.02.  This  would  suggest  the 
pseudospectral  approximation  in  the  vicinity  of  the  wall  was  near  its 
steady  state  value  for  the  above  approach. 
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Examination  of  the  governing  equation  suggests  that  the  first  and 
second  derivatives  and  the  evaluation  of  the  diffusivity  terms  control 
the  solution.  Therefore,  the  validity  of  the  approximations  for 
derivatives  was  examined.  Due  to  the  presence  of  intial  discontinuities 
at  both  boundaries,  the  derivative  approximations  oscillated  about  the 
correct  derivative  values.  Oscillations  were  most  severe  at  the 
boundaries  and  attenuated  as  the  central  region  was  approached.  Since 
the  derivative  values  initially  should  have  been  zero  in  the  core  region 
both  positive  and  negative  deviations  about  the  correct  solution  were 
introduced;  thus  the  digital  filter  that  was  used  for  the  Graetz  problem 
was  applied  to  both  the  first  and  second  derivatives.  This  approach 
produced  reasonable  results  in  comparison  to  the  finite  difference 
values . 

Figures  3  and  4  in  addition  to  Tables  3  and  4  display  the  results 
obtained.  These  results  are  for  a  second-order  Adams-Bashforth 
predictor  with  step  sizes  of  1,  2.5  and  10  seconds  for  33,  25  and  17 
term  representations  respectively.  The  derivatives  were  filtered  for 
the  first  500  seconds,  in  addition  a  "cut-off"  filter  was  used. 
Typically,  due  to  the  oscillation  of  the  derivative  values,  the  solution 
also  oscillated  to  a  lesser  degree,  most  noticeably  when  the  true 
solution  was  constant.  Therefore  when  the  collocation  approximation  for 
mole  fraction  exceeded  one,  that  point  along  with  the  succeeding  points 
were  set  equal  to  one,  thus  "cutting-off "  oscillatory  behavior.  The 
finite  difference  representation  used  nodes  spaced  0.005  cm  apart  and 
symmetry  about  x=0  with  the  integration  being  carried  out  by  a  second- 
order  Adams-Bashforth  method  with  time  increments  of  0.5  seconds. 

38 


DISTANCE  FROM  SLAB  CENTER  (an) 


Figure  3. 


Ethanol  concentration  Profiles  within  a  Slab  Obtained  by  a 
Finite  Difference  Representation  and  Pseudospectral 
Representations  of  17,  25,  and  33  Terms,  for  t-1000  seconds 
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Figure  4.  Ethanol  concentration  Profiles  within  a  Slab  Obtained  by  a 
Finite  Difference  Representation  and  Pseudospectral 
Representations  of  17.  25,  and  33  Terms,  for  t=5000  seconds. 
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Table  3. 

Comparison  of  ethanol  mole  fractions  at  various 
a  slab  at  a  time  of  1000  seconds  for  finite  d: 
and  pseudospectral  methods. 

points  with 
Lfference  (Fl 

0.9239 

x(cm) 
0.7071 

0.3827 

0 

*FD 

0.5357 

0.9281 

0.9984 

1 

%  Error 

0.5982 
11.67 

0 . 8958 
3.480 

0.9910 
0.7412 

1 
0 

*24 

%   Error 

0 . 7094 
32.42 

0.9355 
0.7973 

0.9991 
0.07011 

1 
0 

*32 

%  Error 

0.6157 
14.93 

0.9295 
0.1509 

0.9979 
0 . 05008 

1 

0 
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Table  4. 

Comparison  of  ethanol  mole  fracti 
slab  at  a  time  of  5000  seconds 
pseudospectral  methods . 

ons  at  various  points  with 
for  finite  difference  (FD) 

x(cm) 

0.9239 

0.7071 

0.3827 

0 

*FD 

0.1954 

0.7188 

0.9116 

0.9695 

♦ia 

0.5066 

0.7844 

0.9210 

0.9665 

*  Error 

159.3 

9.126 

1.031 

0 . 3094 

*24 

0.3316 

0.7627 

0.9295 

0.9770 

%   Error 

69.70 

6.107 

1.964 

0.7736 

*32 

0.1980 

0.7233 

0.9151 

0.9708 

*  Error 

1.331 

0.6260 

0.3839 

0.1341 
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Figure  3  displays  the  pseudospectral  and  finite  difference 
approximations  for  t=1000  seconds.  Table  3  compares  error  between 
finite  difference  and  pseudospectral  methods.  It  can  be  seen  that  the 
33  term  series  representation  results  in  the  best  approximation  to  the 
solution.  In  addition  Table  3  substantiates  this  veiw  point  by 
comparison  of  relative  percent  error.  Figure  4  and  Table  4  present  the 
results  for  t=5000  seconds.  Again  it  is  seen  that  the  33  term 
representation  accurately  approximates  the  finite  difference  solution. 

It  is  apparent  from  the  deviations  among  the  approximate  solutions 
that  the  17  and  25  term  expansions  cannot  adequately  represent  the 
initial  solution.  Although,  no  effort  was  made  to  optimize  the 
filtering,  it  is  likely  that  the  derivative  values  require  far  less 
filtering  than  was  applied. 

Diffusion  and  Reaction  in  a  Spherical  Pellet 
The  problem  considered  was  the  diffusion  and  concurrent  reaction  of 
a  substrate  within  a  spherical  catalysis  pellet  with  nonlinear 
Michaelis-Menton  kinetics.  It  is  assumed  that  the  pellet  is  homogeneous 
with  diffusion  and  reaction  symmetric  such  that  only  radial  variation 
need  be  considered.  In  addition,  diffusivity  and  density  of  the  system 
are  assumed  to  be  constant.  The  governing  equation  is: 

as   n  t1     9  /  2  3S.,    VmS  ,   , 

3T  =  D  (-2  3r-  (r  iF))  '  FTs"  <59> 

r  m 

S  is  substrate  concentration,   D  is  the  diffusivity  of  the  substrate 

within  the  pellet,  r  is  radial  position,  V  the  maximum  rate  at  which 

m 
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the  reaction  can  take  place  with  respect  to  the  physical  conditions 

(i.e.    temperature)   and   K    a   kinetic  constant.   The  following 

m 

dimensionless  varibles  are  introduced  into  the  governing  equation 

c  -  S_  .    R  .  £_  .      _  tD_ 

C  "  SK  '    R  "  r   '    T  ~   2 
b         s        r 

s 

S.   is  the  bulk  concentration  of  substrate  and  r  is  the  radial  distance 
o  s 

to  the  pellet  surface.   Which  results  in 


ac  _  a^c   2  3c   _ac 

3r  -  9R2  +  R  8t   "  /8+C 


(60) 


where 


V  r2  K 

a  =  -S-2  ;      8   =  -"• 

b  b 

The  initial  and  boundary  conditions  considered  are: 

C(r,0)  =0;   0  <  r  <  1  (61) 

C(l,t)  -  lj    t  >   0  (62) 

3C(0,t) 

-ir-1   =  0;   t  >  0  (63) 

Due  to  the  nonlinear  nature  of  this  problem,  no  analytic  solution  is 
known  to  exist.  Therefore,  a  finite  difference  method  was  executed  for 
comparison.   The  governing  equation  was  integrated  using  a  second  order 

Adams-Bashforth  method.   The  initial  thirty  increments  were  lxlO-5  in 

size  at  which  point  the  time  increment  was  increased  to  lxl0~  .  The  use 
of  small  initial  time  increments  was  to  damp  the  oscillations  present  in 
the  solution  arising  from  the  initial  discontinuity  and  ease  the 
effects  that  filtering  may  have  on  theultimate  solution.   The  values 
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II 


1 

1 

1 

0.9808 

0.1411 

0.2212 

0 . 9239 

-7.843xl0~3 

0.1226 

0.8315 

2.027xl0~3 

0.04517 

0.7071 

-8.452xl0~4 

0.01841 

0 . 5556 

4.737xl0~4 

6.421x10 

0.3827 

-3.270X10"4 

2.105x10 

0.1951 

2.657xl0~4 

7.139x10 

0 

-1.563xl0~4 

3.347x10 

-3 


-4 


-4 


1 
0.1740 

0.09500 

0.02741 

8.861x10 

2.173xl0~ 

4.382x10 

1.155x10" 

2.265x10" 


-3 


-4 


Table  6.  A  comparison  of  concentration  profiles  within  a  spherical 
pellet  after  seven  time  increments  for  no  filtering  and 
filtering  schemes  I  and  II. 
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chosen  for  a  and  /9  were  7.5  and  1,  respectively;  these  values  indicate 
that  for  a  region  in  which  the  substrate  concentration  is  high  (i.e. 
C — >1),  the  kinetics  are  nonlinear  and  for  the  case  where  the 
concentration  is  low,  first-order  kinetics  hold.  The  filters  are  the 
same  as  those  used  with  the  Graetz  problem  with  the  exception  of 
weighting  on  the  central  point,  which  was  increased  by  one  for  both 
filters.  Two  separate  filtering  approaches  were  used,  the  first  (I) 
filtered  the  first  seven  increments,  the  second  (II)  filtered  at 
increments  seven  and  fourteen  until  no  negative  concentrations 
existed.  It  can  be  seen  from  Table  5  that  significant  differences  can 
occur  from  the  use  of  different  filtering  schemes,  most  notably  in  the 
central  region  where  the  solution  approaches  zero.  The  question  that 
arises,  however,  is  the  extent  to  which  filtering  affects  the  numeric 
values.  Table  6  presents  concentration  values  with  no  filtering  and 
filtering  schemes  I  and  II  at  increment  seven.  It  can  be  seen  there  are 
noticeable  differences  between  the  results  of  the  two  methods.  Note  the 
oscillations  present  without  filtering,  these  oscillations  occur  due  to 
the  presence  of  the  discontinuity  at  the  boundary  and  increase  at  each 
time  step  leading  ultimately  to  numeric  instability. 

External  Separated  Flows 
One  of  the  original  goals  of  this  research  was  the  application  of 
pseudospectral  methods  to  external  separated  flows,  notably  flow  about  a 
pitched  flat  plate.  Although  this  problem  is  extremely  difficult,  it  is 
nevertheless  instructive.  Therefore,  consider  the  two  dimensional 
incompressible  Navier-Stokes  equation  for  flow  past  a  pitched  flat  plate 
expressed  in  terms  of  the  stream  function 
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it**  +  %  h  -  *x  h™  =  ^  (64) 

The  transient  response  of  the  flow  to  a  periodic  forcing  function  is  the 
ultimate  goal,  but  the  first  step  is  to  consider  the  steady  state 
solution.  A  Fourier  spectral  solution  for  the  steady  state  form  of 
equation  (64)  was  presented  by  Mei  and  Plotkin  [1984]  for  confined 
laminar  wakes.  The  use  of  the  Fourier  spectral  method  by  Mei  and 
Plotkin  was  made  possible  by  imposition  of  periodic  boundary  conditions. 
As  a  preliminary  step,  the  Falkner-Skan  equation  was  used  to 
estimate  the  accuracy  of  the  pseudospectral  approximation  for  separated 
flow  present  on  a  divergent  wall.  Based  upon  a  fourth  order  Runge-Kutta 
solution  of  the  Falkner-Skan  equation  {ft   =  -0.125) 

f"  +  ff"  +  /8(l-f'2)  =  0  (65) 

with  boundary  conditions 

f  =  f '  =  0;  rj  =  0 

f ■  -  l;  r?  =  »  (66) 

a  Chebyshev  interpolation  function  was  evaluated  for  11,  21,  and  41  term 
representations.  The  evaluation  was  based  on  a  comparison  of  the  first 
and  second  derivatives  found  using  a  pseudospectral  evaluation  and  those 
obtained  by  the  Runge-Kutta  method.  The  coefficients  for  the 
pseudospectral  approximation  are  obtained  by  Fourier  transformation. 
Three  separate  observations  were  made  from  these  evaluations;  1)  In  one 
trial  the  data  used  for  the  value  of  f(rj)  at  each  collocation  point  was 
rounded  to  three  significant  figures  with  N=20.  This  resulted  in  large 
oscillations  of  the  first  derivative,  as  noted  by  Osher  [1984]  and  Majda 
et  al.  [1978].  2)  Increasing  the  number  of  collocation  points  from  11  to 
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21  to  41  did  not  significantly  affect  the  accuracy  of  the  resulting 
approximations.  3)  The  accuracy  of  a  derivative  decreased  approximately 
one  order  of  magnitude  for  each  increase  in  the  order  of  differentation. 
It  should  be  noted  however,  that  the  decrease  in  accuracy  of  succeeding 
derivatives  may  be  due  to  slight  errors  introduced  by  the  use  of  linear 
Interpolation  to  evaluate  point  values  for  the  pseudospectral 
approximation.  A  compilation  of  these  results  is  presented  in  Tables  7 
and  8. 

Based  upon  the  above,  it  was  decided  to  use  a  steady  state  stream 
function-vorticity  representation  for  the  flow. 

9x  3y   3y  3x  "  U\   2  +   .  2}  (67' 

3x    3y 

ax2  ay 

with  the  boundary  conditions 


° =  -  \H +  H]  (68) 


*(x,0)   =  Vrx(x,0)    =  *  (xt0)   =0         -a  <  x  <  a  (69) 

*{x.tm)   =  *(*»,y)   =  ^free  streaB  (70) 

(j(±»,y)    =  u(x,±a>)    =  0  (71) 

2 
u(x,0)   =-|^~|  -a<x<a  (72) 

ay 

Then  the  stream  function  and  vorticity  can  be  approximated  by: 
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N 
*(x,y  )  =  Z  an(x)Tn(y  )  (73) 

n=0 

N 
o(x,yi)  =  Z  an(x)Tn(yj)  (74) 

This  representation  was  chosen  to  make  use  of  the  increased  accuracy  of 
pseudospectral  approximations  along  the  axis  in  which  the  greatest 
amount  of  change  will  occur;  Gottlieb  et  al.  [1984a]  note  that  a  finite 
difference  representation  is  typically  used  along  the  vertical  axis  and 
a  pseudospectral  representation  along  the  horizontal  axis.  The  solution 
is  split  into  an  upper  and  lower  region  for  above  and  below  the  plate 
and  having  the  numerical  boundary  conditions 

(J  (x,0)  =  <J  (x,0)  a  <  x  <  -a  (75) 

\frU(x,0)  =  \|f  (x,0)  a  <  x  <  -a  (76) 

where  the  superscript  u  and  1   refer  to  the  upper  and  lower  regions 
respectively. 

The  following  computational  scheme  was  then  used  to  evaluate  the 
flow.  The  mesh  nodes  for  the  stream  function  were  initialized  using  an 
equation  given  by  Lamb  [1945]  that  had  been  corrected  for  the  angle  of 
attack  and  transformed  to  Cartesian  coordinates  (Appendix  A).  The 
vorticity  was  initialized  at  zero  everywhere  except  on  the  plate  where  a 
first  order  finite  difference  evaluation  was  performed  based  on  methods 
found  in  Roache  [1972] 

U(x,0)  =  w  %  +  0(An)  (77) 
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where  An  is  the  distance  from  the  wall  (w)  to  the  mesh  node  normal  to 
the  wall  (w+1).  The  following  steps  were  then  followed  in  obtaining  a 
solution: 

1)  Determine   the  pseudospectral  coefficients  to  the  Chebyshev 
polynomials  at  each  (x,y.)  combination. 

2)  Starting  at  the  leading  edge  of  the  flow  field  evaluate  the 
derivatives  with  respect  to  y  for  that  column. 

3)  Equation   (67)  is  solved  with  respect  to  u  for  a  finite  difference 
representation  of  the  x  derivative. 

4)  Equation   (68)  is  solved  with  respect  to  ^  for  a  finite  difference 
representation  of  the  x  derivative. 

5)  Steps  3)  and  4)  are  repeated  for  each  value  of  u>  and  ^   in  the 
column. 

6)  Move  to  the  next  column  and  repeat  steps  2)  through  5). 

7)  Repeat  for  the  second  region. 

8)  Satisfy  the  numerical  boundary  conditions. 

9)  Continue  steps  1)  through  8)  until  convergence  is  achieved. 

In  order  to  estimate  the  accuracy  of  the  above  algorithm  a  central 
finite  difference  algorithm  for  the  stream  function  representation  was 
developed.  The  successive  under  relaxation  algorithm  was  based  on  the 
same  boundary  conditions  as  given  by  equations  (69)-(70). 

Both  the  finite  difference  and  pseudospectral  algorithms  were 
explosively  unstable.  The  fact  that  both  algorithms  were  unstable 
suggests  that  improper  boundary  conditions  are  being  applied.  Either 
boundary  condition  for  the  derivative  of  the  stream  function  on  the 
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plate  may  cause  problems.  These  boundary  conditions  have  to  be 
satisfied  since  each  represents  a  velocity  component  on  the  plate  and 
the  Newtonian  no-slip  condition  requires  that  velocity  be  zero  at  the 
surface.  Evaluating  either  condition  with  a  forward  finite  difference 
approximation  implies  that  the  entire  region  around  the  plate  has  \|f  =  0. 
As  an  aside,  insufficient  numerical  boundary  conditions  were  applied  for 
the  pseudospectral  approximation.  The  first  and  second  derivatives  for 
both  stream  function  and  vorticity  should  have  been  equated  across  the 
boundary . 
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Chapter  IV.   CONCLUSIONS 

The  objective  of  this  work  was  to  assess  the  suitability  and 
computational  efficiency  of  the  application  of  spectral  methods  to 
problems  in  transport  phenomena.  This  was  accomplished  by  the 
evaluation  of  representative  problems  from  each  major  area  of  transport 
phenomena;  momentum,  mass  and  heat  transfer. 

It  was  found  that  Chebyshev  pseudospectral  representations  can 
accurately  evaluate  derivatives  of  smooth  functions  with  few  expansion 
terms.  The  reduction  in  the  number  of  nodes  at  which  a  solution  was 
sought,  along  with  the  use  of  matrices  to  evaluate  derivatives  resulted 
in  significantly  faster  solution  times  than  obtained  with  finite 
difference  representations  of  corresponding  accuracy.  The  presence  of  a 
discontinuity  in  the  solution  domain  was  found  to  be  a  serious  problem 
for  a  number  of  initial  and  boundary  value  problems.  A  discontinuity 
affects  the  entire  solution  domain  due  to  the  global  nature  of 
derivative  evaluation.  The  effects  of  this  problem  can  be  removed  by 
the  use  of  digital  filtering;  such  filters  were  found  to  be  quite 
satisfactory  in  the  removal  of  "noise"  from  the  solution  arising  from  a 
discontinuity.  Therefore,  pseudospectral  methods  should  be  applicable 
to  many  if  not  most  problems  of  interest  in  transport  phenomena. 

There  is  need  to  develop  criteria  for  the  application  of  filters 
to  pseudospectral  solutions,  namely  the  type  of  filter  that  should  be 
applied  and  amount  of  filtering  necessary  to  arrive  at  accurate  results 
for  a  given  type  of  problem.  Additionally,  filters  that  retain  the 
"infinite"  order  accuracy  of  spectral  methods  are  desirable.  Finally, 
for   problems   in  which  a  discontinuity  is  present  due  to  initial 
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conditions,  other  methods  for  initializing  the  solution  should  be 
investigated  so  that  the  discontinuity  can  be  removed  and  the  accuracy 
of  the  pseudospectral  method  exploited.  An  example  of  alternative 
initialization  would  be  the  use  of  a  high  order  finite  difference  method 
to  remove  the  discontinuity  and  subsequent  use  of  an  accurate 
interpolation  technique  to  provide  solution  values  for  the 
pseudospectral  points. 
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APPENDIX  A 
STREAMLINE  EXPRESSION  FOR  A  PITCHED  FLAT  PLATE 
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Lamb  [1945]  gives  the  following  representation  for  the  stream 
function  for  flow  past  a  flat  plate  pitched  at  angle  of  45° 

t  = sinh£(cosn  -  sinn)  (Al) 

whose  edges  are  located  at  x  =  ±  c  and  where  q0  is  the  free  stream 

velocity.   This  equation  was  derived  from  the  relation  describing  fluid 
motion  relative  to  an  elliptic  cylinder 

^  =  -  c  J  v  sinh^cosn  -  u  sinh£sin£J  (A2) 

where  u  and  v  are  the  x  and  y  components  of  velocity  with  respect  to  the 
plate.   Then  \^  can  be  represented  as 

i/   =  -c  q0cos8sinh?|tanecosn-sinrj|  (A3) 

where  ©  is  the  angle  of  attack  for  the  plate.   Cartesian  coordinates  and 

elliptic  coordinates  are  related  in  the  following  fashion 

x  =  ccosh(£)cos()7)  (A4) 

y  =  csinh(£)sin(n)  (A5) 

The  following  relation 

2         2 
x         y 

i — I — s — r  =  l  (A6) 

c  cos  r)   c  sin  n 
along  with  the  definition  of  a  point  on  a  hyperbola  and  the  relation  of 
sinh  to  cosh  define  the  following  relations  for  the  stream  function  in 
Cartesian  coordinates 


*  =  -q0cose{  tane\A2  -  |(  >/(c+x)2+y2'  -  \/(c-x)2+y2'  )2  -y  }   (A7) 


for  y>0,  and 


*  =  q0cose{  taneyk2  -  \(   v^c+x)2+y2'  -  \Xc-x)2+y2' )2  -y  } 


(A8) 


for  y<0. 
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APPENDIX  B 
COMPUTER  PROGRAMS 
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The  programs  in  this  Appendix  represent  each  section  of  Chapter 
III.  Two  programs  have  been  used  for  the  evaluation  of  the  derivative 
matrix  defined  by  equation  (32)  given  N,  where  there  are  N+l  terms  in 
the  series  expansion  (0  to  N).  All  programs,  with  the  exception  of 
one,  have  been  written  in  Basic  due  to  the  availability  of  interpative 
and  compiler  versions  of  this  language  for  microcomputers.  The  one 
exception  was  written  in  Fortran.  The  programs  should  be  self 
explanatory  with  the  remarks  provided  in  each. 
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C  THE  FOLLOWING  PROGRAM  EVALUATES  MATRICES  FOR  THE 

C  DETERMINATION  OF  DERIVATIVE  COEFFICENTS. 

REAL  A  (  33 , 33 )  ,  E  ( 33 , 33 ) , C ( 33 , 33 ) , D ( 33 , 33 ) 
DATA  A, E, C, D/ 10 89*0. , 10  89*0. , 1089*0.  ,  1089*0. / 
DATA  L, M, N,  IA,  IE,  IC/33, 33, 33, 33, 33, 33/ 
OFEN(l, FILE='C1.  DAT'  ) 
QPEN(2,FILE='C2.DAT' ) 

C      OPEN (3,  FILE='C3.  DAT' ) 

C      0PEN(4,FILE='C4.  DAT' 
WEITEi*, 753) 

753  FORMAT ( '  N? ' ) 
READ U, 754 )N 

754  FORMAT<I2) 
INDEX =2 
AINDEX=N-2 
EINDEX=N-2 
DO  10  1=1, N 

TI  =  I 

DO  20  J=INDEX,N, 2 

A (I, J)=  2.  *(N-1. -B INDEX) 
EINDEX=BINDEX  -  2. 
IF  (J.EQ.N-1)  GOTO  16 
GOTO  20 
16  IF  (2*(I/2).NE. I)  GOTO  15 

20  CONTINUE 
15  INDEX=INDEX+1 

IF  ( INDEX. GT.N)  GOTO  26 
EINDEX=AINDEX-TI 
10  CONTINUE 

26  CALL  VMULFF(A, A, L, M, N,  IA,  IE, B,  IC, IER) 
CALL  VMULFF(A,B,L,  M,N,IA,IB,C,IC,IER) 
CALL  VMULFF(B,B, L, M, M, I A, IE, D, IC, IER ) 
DO  30  1=1, N 

A(l, I)=A(1, I)/2, 

E(l, I)=B(1, I)/2. 

C(l, I)=C(1, I)/2. 

30      D(l,  I)=D(1,  D/2. 

DO  90  1=1, N 

DO  90  J  =  1,N 

WRITEd,  1000)A(I,  J) 
WRITE (2, 1000 )B( I,  J) 
C  WRITE (3, 1000 )C ( I,  J) 

C  WRITE (4, 1000 )D( I, J) 

90  CONTINUE 
1000  F0RMAT(F15. 2) 
STOP 
END 


61 


i  '  THE  FOLLOWING  PROGRAM  DETERMINES  THE  MATRIX  FOR  THE 
2  •'  DIRECT  EVALUATION  OF  DERIVATIVES  FROM  SOLUTION  VALUES 
10  CLS 

15  FI#=3. 141592654* 

16  N  =  3 2 

20  DIM  01(33, 33) , 02(33, 33) , Tl (33, 33) , T2(33,  33  5 

21  DIM  TT( 33, 33), TO (33,  33) 
30  OPEN" I", 1, "CI. DAT" 

40  OPEN" I", 2, "C2. DAT" 
50  FOR  1=0  TO  N 
60      FOR  J  =  0  TO  N 
70  INPUT#1,C1( I, J) 

80  IMPUT#2, C2(I, J) 

90      NEXT  J 
100  NEXT  I 
110  FOR  1=0  TO  N 

120      IF  1=0  OR  I=N  THEN  01=2  ELSE  Cl=l 
130      FOR  J=0  TO  N 

140  IF  J=0  OR  J=N  THEN  C2=2  ELSE  02=1 

150  TT(I,J)=  2*COS( I*J*PI#/N)/N/C1/C2 

160  TC( J, I)=COS(I*J*PI#/N) 

170      NEXT  J 
130  NEXT  I 
l'?0  FOR  1  =  0  TO  N 
200      FOR  J=0  TO  N 
210  Tl( I,J)=0: T2( I, J)=0 

220  FOR  K=0  TO  N 

230  Till,  J)=T1( i, J)     +  Cl( I,K)*TT(K, J) 

240  T2(I, J)=T2( I, J)  +  C2(I,K)#TT(K, J) 

250  NEXT  K 

260      NEXT  J 
270  NEXT  I 
230  FOR  1=0  TO  N 
*290      FOR  J  =  0  TO  N 
300  OKI,  J)=0:C2(  I,  J)=0 

310  FOR  K=0  TO  N 

320  CKI,  J)=C1(I,  J)  +  TC(I,K)*T1(K, J) 

330  C2(I,J)=C2(I,J)  +  TC(I, K)*T2(K, J) 

340  NEXT  K 

350      NEXT  J 
360  NEXT  I 

370  0?EN"0",3, "G1N32. DAT" 
380  OPEN"0", 4, "G2N32. DAT" 
390  FOR  1=0  TO  N 
400      FOR  J=0  TO  N 
410  ?RINT#3,C1(I, J) 

420      NEXT  -J 
430  NEXT  I 
440  FOR  1=0  TO  N 
450      FOR  J=0  TO  N 
460  FRINT#4, 02(1, J ) 

470      NEXT  J 
43  0  NEXT  I 


490  CLOSE 

1  '  The  following  program  calculates  the  start-up  velocity 

2  '  profile  in  a  circular  tube  using  Chebyshev  pseudo- 

3  '  spectrai  method. 
10  CLS 

20  DIM  V(25),G1(25,  25),G2(25,  25) , VP ( 25 ) , V?P < 25 ) ,FGLD(25) 

30  N=24:PI=3. 1415926* 

40  0PEN"I",1,  "G1N24.DAT" 

50  QFEN"I", 2, "G2N24.DAT" 

60  FOE  1=0  TO  N 

70      FOE  J=0  TO  N 

80  IN?UT#1, Gl( I, J) 

90  INPUT#2.G2(I, T) 

100      HEXT  J 

110  NEXT  I 

120  FOR  JFK=1  TO  800 

130      GOSUE  300 

14  0      FOR  1=1  TO  N/2 

150      R=C0S(I*PI/N): IF  R<. 00001  THEN  R=. 00001 

160  V(I)=V(I)  +  .000125* (4  +  VP(I)/R  +  VPP(D) 

170      NEXT  I 

180      FOR  I=N/2  TO  N 

190  V(I)=  V(N-I) 

20  0      NEXT  I 

210      FOR  1=0  TO  N/2 

220  PRINT  V( I) 

230      NEXT  I 

240  NEXT  JFK 

250  LPRINT  "N=24    UNSTEADY  LAMINAR  FLOW  IN  A  CIRCULAR  TUEE" 

260  FOR  1=0  TO  N/2 

270      LPRINT  I,COS(I*PI/N) , V(I) 

2  80  NEXT  I 

290  STOP 

300  FOR  1=1  TO  N/2 

310      VP(I)=0: VPP(I)=0 

320      FOR  u=l  TO  N 

330  VP(I)=VF(I)  +  G1(I,J)#V(J) 

340  VPP( I )=VFP ( I )  +  G2( I, J)*V( J) 

350      NEXT  J 

360  NEXT  I 

370  RETURN 
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i    The  following  program  evaluates  the  transient 

2  temperature  profile  of  a  wire  electrically  heated. 

3  The  Chebyshev  pseudcspectrai  method  is  used  in  this 

4  '  program. 

10  CLS  .   . 

20  DEFINT  I 

30  DEFINT  J 

4  0  N  =  24 

50  DIM  7(60) , TFP(60 ) , A (60 >, AFP (60) , OLD (60) 

60  DIM  G2(25, 25) 

70  0FEN"I\1,  "G2N24.DA7" 

SO  FOR  1=0  70  N 

90      FOH  J=0  70  N 

100  INPUT#1,G2(I, J) 

110      NEXT  J 

120  NEXT  I 

130  H=.0001 

140  PI =3. 1415926* 

150  ALFHA=. 1 

160  EETA=8.5317E-11 

170  7M=i.95555E+12 

180  70=281 

190  FOH  1=0  70  N 

200     7(1) =70 

210     0LD( I )=167. 28 

220  NEX7  I 

230  LFRIII7  TIME* 

240  FOR  JFK=1  70  400 

250  GOSUE  410 

260      FOE  1=1  70  N/2 

270  7EMF=  ALPHA*  7PP(  I)  +  EE7A*(7M  -  7(1) -•4) 

280  7(1)=  7(1)  +  H*(3*7EMF  -  0LD(I))/2 

290  0LD(I)=7EMF 

300      NEX7  I 

310      FOR  1=1  70  N/2 

320  7(N-I)=7(I) 

330      NEXT  I 

340  NEXT  JFK 

350  LFRINT  TIME* 

360  OF EN "O", 2, "NIRE5.DAT" 

370  FOH  1=0  70  N/2 

33  0  FRINT#2,C0S(I*PI/N)f T( I), A*, C0S( ( I +N/2+1 ) *PI /N ) , 
Td  +  N/2+1) 

390  NEX7  I 

400  370P 

410  FOR  1=1  70  N/2 

420      7FF(I)=0 

430      FOR  1=0  70  N 

440  7FP( I )=7?P (  I  )  +  7( J) *G2( I, J) 

450      NEX7  J 

460  NEXT  I 

470  HE7URN 
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1  The  following  program  evaluates  the  classical  Graetz 

2  '  problem  for  ?e=l  and  He=l. 
10  DEFDBL  C 

20  C  L  S 

30    DIM    C(25)  ,C?(25),  CP?<25>,  COLD  (25),  GK25,  25) 

40    DIM    C10LD(25) , G2<25,  25) , CM (25) 

50  PI  =  3. 141592654# 

60  OPEN "0", 3, "GEAETZ.DAT" 

70  N=16 

80  ALPHA=7.5 

90  EETA=1 

10  0  FE=1 

110  H=. 00001 

120  FOE  1=1  TO  N-l: C(I)=1: NEXT  I 

130  OPEN" I" , 1, "G1N16. DAT 

140  OPEN" I", 2, "G2N16. DAT 

150  FOR  1=0  TO  N 

160      FOR  J=0  TO  N 

170  INPUT#i,GKI,  J) 

180      NEXT  J 

190  NEXT  I 

200  FOR  1=0  TO  N 

210      FOR  J=0  TO  N 

220  INPUT#2, G2< I, J)   ■ 

230      NEXT  J 

240  NEXT  I 

250  '***##**##***########*##*###*#**#***#* #****#***#*#**#* 

260  C(0)=0: C(N)=0 

270  FOR  JFK=1  TO  7501 

280      GOSUE  650 

290      FOR  1=1  TO  N/2 

300        -  R=COS( I*PI/N) 

310  IF  I=N/2  THEN  R=. 00005 

320  TEMP=  (CPP(I)  +  CP(I)/R)/(1  -  EA2)/PE 

330  C(I)=  C(I)  ♦  H*(3*TEMP  -  C0LD(I))/2 

340  COLD(I)=TEMP 

350      NEXT  I 

360      FOR  1  =  0  TO  N/2:  'PRINT  COS ( I*PI/N) , C ( I ) 

370      C(N-I)=C(I) 

380      CN(I)=C(I) 

390      NEXT  I 

400  'IF  JFK=20  THEN  GOSUE  670 

410  'IF  JFK=30  THEN  GOSUB  710 

420   IF  JFK=2500  OR  JFK=5000  OR  JFK=75O0  THEN  GOSUE  610 

430   IF  JFK>7  THEN  GOTO  580 

440      CN ( N/2  +  1 ) =C (N/2-1  ) 

450      CN(N/2+2)=C (N/2-2) 

460      C(l)=  (11*CN(1)  +  CN(2))/12.25 

470      FOR  1=2  TO  N/2 

480  R=COS ( I*FI/N) 

490  RP1=C03( ( 1+1 ) *PI/N) 

500  EF2=C0S( ( 1+2) *FI/N) 

510  RM1=C0S( (I-1)*PI/N) 


520  RM2=COS( (I-2)*PI/N) 

530  .        Xl=  2*<R  -'HM1)/(RP1  -  R) 

540  X2=(RM2  -  R)/(H  -  RF2) 

550  Cm  =  (X2*CN(I-2)+Xl*CN(I-l.)+10*CN(I  )+2#CN<  1  +  1! 

+  CN(  1  +  2)  )/'(13  +  XI  +  X2) 
560      NEXT  I 

570  FOR   1  =  0   TO   N/2:  PRINT   COS < I *P I /N) , C i I ) , CN ( I  ) 

571  C(N-Ii=C(I) 
575     NEXT  I 

58  0  NEXT  JFK 

59  0  CLOSE 

60  0  STCF 

610  FRINT#3,JFK 

620  FOR  1  =  0  TO  N/2: ?RIMT#3, COS ( I *?  I  /N ) , C ( I ) : NEXT  I 

630  FEINT#3, "  " 

640  RETURN 

650  FOR  1=1  TO  N/2 

660      CF(I)=0 

670      CPP(I)=0 

68  0      FOR  J=0  TO  N 

690  CF(I)=CF(I)  +  Gl( I, J)*C(J> 

700  CFF( I)=CFF< I )  +  G2(I, J)*C(J) 

710      NEXT  J 

720  NEXT  I 

730  RETURN 

740  FOR  1=0  TO  N 

750      C10LD< I)=COLD( I ) 

760  NEXT  I 

770  RETURN 

780  FOR  1=0  TO  N 

790      COLD (I )=C10LD( I ) 

800  NEXT  I 

810  H=. 00005 

820  RETURN 
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1  ■'  THE  FOLLOWING  PROGHAM  EVALUATES  THE  DIFFUSIVITY  IN  A 

2  '  POROUS  PLANAR  SLA3.   THE  FIRST  AND  SECOND  DERIVATIVE 

3  '  ARE  FILTERED  WITH  A  DIGITAL  FILTER. 
10  DEFINT  I 

20  DEFINT  J 

30  CLS 

40  DIM  C(33), CP(33) ,CFP(33), C0LD(33) ,01(33,33) 

50  DIM  C10LD(33), 02(33, 33) , CN(33) 

60  PI=3. i4I592654# 

70  N=24 

SO  A=i. 443935 

90  B=-7. 373196 

100  CA=16. 15581 

110  D=-7.9731B4 

120  H=2 

130  OPEN" I " , 1, "G1N24. DAT 

140  OPEN" I", 2, "G2H24.DAT 

150  FOR  1=0  TO  N 

160      FOR  J  =  0  TO  M 

170  INPUT#1,G1(I, J) 

180      NEXT  J 

190  NEXT  I 

200  FOR  1=0  TO  N 

210      FOR  J  =  0  TO  N 

220  IHPUT#2, G2( I, J) 

230      NEXT  J 

240  NEXT  I 

250  '*********#****#***#****##*#***********#*******###*•*** 

260  FOR  1=1  TO  N-l:C( I)=l: NEXT  I 

270  FOR  JFK  =  1  TO  500 

280      GOSUE  670 

290      FOR  1=1  TO  N/2 

300  DAE=  (A  +  E*C(I)  +  CA*C(I)A2  +  D«C < I ) A3 J *. 00001 

310  DAEX=  (E  +  2*CA*C(I)  +  3*D*C ( I ) A2 ) * . 00001 

320  TEMP=  DAEX*CP(I)A2  +  DAE*CPP(I) 

330  C(I)=  C(I)  +  H*(3*TEMF  -  COLD(I))/2 

340  COLD(I)=TEM? 

350      NEXT  I 

360  FOR  1=0  TO  N/2: 'PRINT  COS ( X*PI/N) , C ( X ) 

370  IF  C(I)>1  THEN  CHK  =  2 

380  IF  CHK  =  2  THEN  C(I)=1 

390  PRINT  COS( I*PI/N) , C( I  ) 

40  0  C(N-I)=C(I) 

410  CN(I)=C(I) 

420  NEXT  I 

430  CHK  =  1 

440  'IF  JFK=98  THEM  GOSUE  670 

450  'IF  JFK=100  THEN  GOSUE  710 

460  IF  20*CINT(  JFK/20  )0  JFK  THEN  GOTO  620 

470  IF  JFK>0  GOTO  620 

48  0  CH( N/2+1 )=C( N/2- i) 

490  CN( N/2+2) =C (N/2-2) 

500  C(l)=  (.25  +  11*CM(1!  +  CN(2))/12.25 


510  FDR  1=2  TO  N/2 

520  R=C03(I*PI/N) 

530  RP1=C03( ( I+1)*PI/N) 

540  R?2=C03( (I+2)*PI/N> 

550  RH1=CQ3< i I-1)*FI/N) 

560  RM2=C03( ( I-2)*?I/N> 

570  XI=  2*(R  -  RM1)/(RF1  -  R) 

530  X2=  (HM2  -  R)/(R  -  RP2) 

590  C(I  )=  (X2*CNCI-2)+Xl*CN(I-l)+10*CN(I)+2*CN(I  +  l) 

+  CN(I+2) )/(13  +  XI  +  X2) 

60  0  NEXT  I 

610  FOR  1=0  TO  N/26 

611  PRINT  C03( I*?I/N), C(  I), CH(  I) 

612  C(N-I)=C(I) 

613  NEXT  I 
620  NEXT  JFK 

630  0FEN"0",3, "5LAESH. DAT" 

640  FOR  1=0  TO  N/2: PRINT#3, COS ( I *P I /N ) , C ( I ) : NEXT  I 

650  CLOSE 

660  STOP 

670  FOR  1=1  TO  N/2 

680  C?(I)=0 

690  C?F(I)=0 

700  FOR  J=0  TO  N 

710  CF(I)=C?(I>  +  GKX,  J)*C(J) 

720  C?P(I)=CP?(I)  +  G2(I,J)*C(J> 

730  NEXT  J 

740  'IF  AES(CPP(I)/CFP(1) X. 01  THEN  CPP(I)=0 

750  'IF  AES(CP(  I)/CP(1)  )•;:.  01  THEN  CF(I)=0 

760  NEXT  I 

770  IF  JFK>250  THEN  RETURN 

730  FOR  1=0  TO  N/2: CN ( I ) =CP ( I ) : NEXT  I 

790  CF(1)=  (14*CN(1)  +  CN(2))/15 

3-00  FOR  1  =  2  TO  N/2-2 

810  R=COS( I*FI/N) 

820  RP1=C0S( ( I+1)*PI/N) 

830  RF2=C08( (I+2)*PI/N) 

340  RM1=C0S( (I-1)*PI/N) 

850  RH2=COS( (I-2)*?I/N) 

360  Xl=  2*(R  -  HM1)/(RP1  -  R) 

870  X2=  (HM2  -  R)/(R  -  RP2) 

880  CP(I)=  (X2*CN( I-2)+Xl#CN(I-l)+13*CN(I)+2*CN( 1+1) 

+  CN( 1+2) )/( 16+X1+X2) 

890  NEXT  I 

900  FOR  1=0  TO  N/2: CN( I ) =CP?( I ): NEXT  i 

910  CP?(1)=  <14*CN(1)  +  CN<2))/15 

920  FOR  1=2  TO  N/2-2 

930  R=COS(I*PI/N) 

940  RP1=C0S( (I+1)#PI/N) 

950  RP2=COS ( ( 1+2) *?I/N) 

960  RN1=C0S £ ( I-1)#PI/N) 

970  RM2=COS( (I-2)*?I/N5 
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980  Xi=  2*<R  -  EM1)/(R?1  -  R) 

990  X2=  (RM2  -  R)/(R  -  RP2) 

10  00        CFPm  =  (X2*CN(I-2)+Xl*CN(I-l)+13*CN<I5+2*CN(  1  +  1) 

+  CN( 1+2) )/<16+Xl+X2) 

1010      NEXT  I 
10  20  RETURN 


1  The  following  program  evaluates  the  transient  concen- 

2  '  tration  profile  within  a  spherical  nodule  in  which 

3  '  diffusion  and  reaction  is  occuring  simultaneously. 

4  '  The  reaction  kinetics  considered  are  those  proposed 

5  by  Michaeiis-Menton  for  enzymes. 
10  CLS 

20  DIM  C(25),C?(25), CP?(25), C0LD(25>, G1C25, 25) 

30  DIM  CN(25),G2(25, 25), C1QLD(25) 

40  PI=3. 141592654* 

50  H=16 

60  AL?HA=7.5 

70  EETA=i 

30  H  =  . 00001 

90  OPEN" I" , 1, "G1N16. DAT 

100  OPEN" I" ,2, "G2N16.DAT 

110  FOR  1=0  TO  N 

120      FOR  J=0  TO  N 

130  INFUT*1,G1(X, J) 

140      NEXT  J 

150  NEXT  I 

160  FOR  1=0  TO  N 

170      FOR  J=0  TO  N 

130  INPUT#2,  G2<  I,  .J) 

19  0      NEXT  J 

200  NEXT  I 

220  C(0)=l:C(N)=l 

230  FOR  JFK=1  TO  277 

240  GOSUE  530 

250  FOR  1=1  TO  N/2 

260  R=CQS(I*FI/N) 

270  IF  I=N/2  THEM  R=. 00005 

280  TEMP=CPP(I )+2*CP(I )/R-AL?HA#C( I)/(BETA+C( I) ) 

290  C(I)=  C(I)  +  H*(3*TEMP  -  CQLD(I))/2 

300  -COLD(I)=TEMP 

310  NEXT  I 

320  FOR  1  =  0  TO  N/2:  PRINT  COS ( I *P I/N ) , C  (  I  ) 

330  C(N-Z)*C(Z) 

340  CN(I)=C(I) 

350  NEXT  I 

360  IF  JFK=20  THEM  GOSUE  670 

370   IF  JFK=30  THEN  GOSUE  710 
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330  IF  7*CINT(JFK/7)=JFK  THEN  GOTO  331  ELSE  GOTO  530 

331  IF  JFKM4  GOTO  530 
385  I!;DEX  =  0 

390      CN(N/2+l)=C< N/2-1) 

400      CH(N/2+2)  =C  (N/2-2) 

410      C(l)=  (.25  +  11#CN(1)  +  CN(2))/12.25 

420      FOE  1=2  TO  N/2 

430  R»CQS(I*PI/N) 

440  RP1=C0S( (1+1) *PI/N) 

450  EP2=COS ( ( 1+2) *PI/N) 

460  RM1=C0S( < I-i)*PI/N) 

470  RM2=C0S( <I-2)*PI/N) 

480  Xi  =  2*(R  -  RM1)/(R?1  -  E) 

490  X2=  (RM2  -  E)/(E  -  RF2) 

500  C( I)=(X2*CN( I-2)+Xl*CN(I-l)+10*CN(I >+2*CN( 1+1) 

+  CN< 1+2) i/(13+Xl+X2i 
505  IF  CdKO   THEN  INDEX=1 

510      NEXT  I 
5i5      IF  INDEX=1  THEN  GOTO  320 

520  FOR  1=0  TO  N/2 

521  PRINT   COS( I*PI/N) , C(I),CN(  I) 

522  C(N-I)=C(I) 
5i5     IitAi  I 

530  NEXT  JFK 

540  0?EN!,0",3,  "NODP.DAT" 

550  FOR  1  =  0  TO  N/2: FRINT#3, COS ( I *P I/N  )  , C C I  )  : NEXT  I 

560  CLOSE 

570  STOP 

580  FOR  1=1  TO  N/2 

590      CP(I)=0 

60  0      C??(I)=0 

610      FOR  J=0  TO  N 

620  CP(I)=C?( I)  +  Gl( I, J)*C( J) 

630  CPP( I ) =CPP( I )  +  G2( I, J)*C(J) 

640      NEXT  J 

650  NEXT  I 

660  RETURN 

670  FOR  1=0  TO  N 

680      C10LD( I )=COLD( I) 

690  NEXT  I 

70  0  RETURN 

710  FOR  1=0  TO  N 

720      COLD(I)=C10LD(  I) 

730  NEXT  I 

740  H=.00  01 

750  RETURN 
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1'  THE  FOLLOWING  PROGRAM  ATTEMPTS  TO  CALCULATE  FLOW  PAST  A 

2'  PITCHED  FLAT  PLATE. 

10  CL3 

20  DIM  PC42, 26) , ?X(42, 26) , FY(42, 26) , FXX(42, 26) , ?7Y (42, 26) 

30  DIM  ?XXX(42, 26),PYYY(42, 26 ) , PXXY Y ( 42, 26), PXXY(42, 26) 

40  DIM  PYYT(42) ,?YYYT(42), ?YYX(42, 26) 

50  'OPEN" I", #1, "FPLATE.DAT" 

60  IMFUT  "NUMBER  OF  ITERATIONS "; HUM 

70  'FOR  1=1  TO  42 

80  'FOR  J  =  l  TO  26 

90  'INFUT#1,P(I, J)   | 

10  0  'NEXT  J : NEXT  I 

110  'CLOSEtl 

120  HUM  =  mJM-l 

130  CLS 

140  MIT  =  0:  H  =  . 6666667: VIS  =  .  17 

150  'GOTO  410 

160  '********  INTIALIZE  NODAL-  POINT  ************************ 

170  FOR  1=1  TO  42 

180  FOR  J=l  TO  26 

190  P(I,J)=  -1.5625*1  +  3.125*J  -  3.90625 

200  IF  J>9  THEN  P  (  T,  J  )  =-1.  5625*  I  +3. 125*J  +  3.90625 

210  NEXT  Jl  NEXT  I 

220  '*******  INTIALIZE  POINTS  ON  THE  FLATE  ***************** 

230  J  =  9 

240  FOR  1=11  TO  26 

250  P(I,J)=0:PY(I,J)=0:?X(I,J)=0 

260  NEXT  I 

270  '*****  AVERAGE  THE  VALUES  OF  THE  STREAM  FUNCTION  ON  ft*** 

280  '***********  THE  RIGHT  AND  LEFT  BOUNDEIES  ************** 

290  FOR  1=1  TO  42  STEP  41 

300  FOR  J=3  TO  24 

310  P(  I,  J)=(?(  I,.J-i)  +P(  I,  J)  +?(  I,  J  +  l)  )/3 

320  NEXT  J : NEXT  I 

330  FOR  1=2  TO  10 

340  FOR  j=3  TO  24 

350    P(I, J)=(P(I, J-2)+P(I, J)+P(I, J-1)+?(I, J+1)+P(I, J+2) )/5 

360    P ( 1  +  25, J)  =  (P ( 1  +  25,  J-2)+P ( 1  +  25, J) +P( 1  +  25, J-l) 

+P( 1+25, J+l) +P( 1+25, J+2) )/5 
370  NEXT  J 
38  0  NEXT  I 
390  FOR  1=  16  TO  26 
400  P(Ii-2,7).=  P(I  +  2,7)/3- 
410  ? ( I, 8) =  P( I, 8  )/12 
420  'P(  I, 10)=P(I, 10)/12 
430  'P(I  +  2,11)=  P(I+2,  11) *(29-I )/12 
440  NEXT  I 
450  FOR  1=27  TO  32 
460  P(I,S)=  F(I,3)*(I-26)/12 
470  ?(I,9)=  P( I, ?) *( 1-26) /12 
480  ?(I,7)=  P( I, 10) *( 1-26) /3 
49  0  NEXT  I 
500  FOR  K=5  TO  19 


510  PRINT  USING"###. ###  "; ?(20, 23-K) , Fi2i, 23-K), P<22, 23-K) 
, P(23, 23-K), F (24,  23-K),  P(25, 23- K),P C 26,  23-K > , p ( 27, 23-K ) 
, ? ( 23, 23-K) 

w ■  i.  0         iJtAl    t'. 

530  F03  1=3  TO  40 

540  FOE  7=3  TO  24 

550  '****  CALCULATE  THE  1st,  2nd,  Si  3rd  DERIVATIVES  ******* 

560  I?  J=  9   AND  I>10  AND  K27  THEN  GOTO  590 

570  PX(I,.J)=  (P(I+1,J)  -  P(I-1, J) J/2/H 

580  FY(I,J)=  (P(I,J+1)  -  P(I, J-l) )/2/H 

590  PXX(I,J)=  <P(I+i,J)  -  2*P(I,J)  +  ?(I-1, J) )/HA2 

60  0  PXXXJ I,J)=(P(I+2,J)-2*P(I+l,J)+2*F(I-i,J) 

-?( 1-2, J) >/2/HA3 
610  IF  j=  9  AND  1)10  AND  K27  THEN  GOTO  710 
620  FYYCI,J)=  (PiI.J+1)  -  2*F(I,J)  +  ?(I, J-l) )/HA2 
630   IF  J  =  S  AND  I? 11  AND  K26  THEM  GOTO  670 
640   IF  J  =  i0  AND  I;  11  AND  K26  THEN  GOTO  690 
650   FYYY( I, J)  =  (P (I, J+2)~2*?( I , J  +  l  ) +2*? ( I , J-l) 

-?( I, J-2) )/2/HA3 
660  GOTO  750 

670  PYYY(I,J)=  (  P( I, J+1)-3*P(I, J)+3*P(I, J-l>-P« I, J-2) ) /HA3 
68  0  GOTO  750 

690  PYYYd,  J)=  (  -P(  I,  J-1)+3*P(  I,  J)-3*P(I,  J  +  l)+P(  I,  J+2)  )/H''3 
700  GOTO  750 

710  PYYT(I)  =  (P( I,  J  +  2)  -  2*P(I,J  +  1)  +  P(I,J))/HA2 
720  FYYYT(I)=  (-F(I,J)+  3*P  !  I ,  J+l ) -3*P  (  I,  J  +  2  ) +F  (  I ,  J  +  3  )  ) /H''3 
730  PYY(If J)=(?( I, J-2)  -  2*P(I,J-1)  +  P(I,J))/HA2 

PYYYf I, J)=  (P( I, J)-  3*? (I, J-1)+3*P(I,  J-2)-P( I, J-3)  ) /HA3 

NEXT  J 

NEXT  I 

'***********#**  CALULATE  THE  CROSS  DERIVATIVES  ****** 

FOR  1=3  TO  40 

FOR  J=3  TO  24 

IF  J=  9  AND  I>10  AND  K27  THEN  GOTO  850 

IF  J=10  THEN  GCSUB  1300 

PXXYd,  J)=  (PXX(  I,  J+1)-FXX(I,  J-l)  )/2/H 

PYYX(I,J)=  (PYY(I+1,J)-FYY(I-1,J) J/2/H 

PXXYY( I,  J)=  (PXX( I, J  +  1)-2*PXX( I, J)+?XX( I, J-l)  )/HA2 

NEXT  J 

NEXT  I 

'****»****  ITERATIVE  EQUATION  FOR  PSI  ************** 

FOE  J=3  TO  24 

FOR  1=3  TO  40 

IF  J=  9  AMD  I>10  AND  K27  THEN  GOTO  10  80 

A=(PY(I,  J)*(PYYX(I, J)+PXXX(I, J)  )-PX(I, J)*(PYYY( I, J) 
+PXXY( I, J) ) )*HA4/12/VIS 

B=-HA4*PXXYY(I, J)/6 
IF  J=  8  AND  I>11  AND  IC26  THEM  GOTO  980 
IF  J  =  10  AND  I>11  AND  K26  THEN  1030 

C=-. 0  3333* (?( I-2,J)+?(I,J-2)+P(I+2,J)+F(I,J+2)) 

D=. 3333* (P( I-1,J)+?(I,J-1)+P(I+1,J!+P(I,J+1)) 

GOTO  1070 

A=   A*12/7 


740 

750 

760 

770 

78  0 

790 

800 

810 

8  20 

830 

840 

850 

860 

370 

880 

390 

90  0 

910 

9  20 

930 

940 

950 

960 

970 

980 

?90  E=  B*12/7 

1000  C  =  -;F(I,J-4!  +  Fil-Z'.J)  +  P(I+2,J))/7 

1010.3=   •+*(?(  1-1,  J)  +  P(I  +  1,.J)  +  P(I,,J-I)  +?(  i,  j-3)  )/? 

-  b  *  P  i  I ,  J  -  2  )  /'  7 
10  20  GOTO  1070 

1030  A  =   A*12/7 

1040  B  =   B*12/7 

1050  C=  -<F(I,J"  +  4)  +  F(I-2,J)  +  P(I+2,J))/7 

1060  D=   4#(P(I-lfJ)  +  P(I+1,J)  f  Pd.J  +  l)  +P(I,J+3))/7 

-  6*P(I,J+2)/7 

1070  P( Ir J)=(A+B+C+D+90*P( I, J) )/100 
10  8  0  NEXT  I 
1090  NEXT  J 
1100  J=2 

1110  FOP  1=2  TO  41 

1120  P(I, J)=(F(I, J-l)+P( I, J+l) )/2 
1130  NEXT  I 
'1140  J=J>23 

.1150  IF  J>26  THEN  1160  ELSE  1110 
1160  NIT=NIT+1 

1170  IF  NIT>NUM  THEN  1180  ELSE  530 
1180  FOR  K  =  5  TO  19 
1190  PRINT  USING"####. ###  " ; P < 20, 23-K ) , P ( 21, 23-K ) , P ( 22, 23-K ) 

, ?(23, 23-K)  ,  F(24,  23-K) ,P(25, 23-K) ,P(26, 23-K) , P ( 27,  23-K) 

, P (28, 23-K) 
1200  'LFRINT  USING  "##.###  " ; F ( 5, 17-K) , P ( 6, 17-K ) , P (  7,17-K) 

,?(8, 17-K),P(9, 17-K), P( 10, 17-K), P ( 11, 17-K ) , P ( 12, 17-K) 
,P(13, 17-K),F(14, 17-K),? (15, 17-K) 
1210  NEXT  K 

1220  OFEN"0", #1, "FPLATE. DAT" 
1230  FOR  1=1  TO  42 
1240  FOR  J=l  TO  26 
1250  FRINT#i,P(I, J) 
1260  NEXT  J 
1270  NEXT  I 
1280  CLOSE 
1290  END 

1300  FOR  L=ll  TO  26 
1310  PYY(L,  9)=PYYT(L) 
1320  ?YYY(L,  9)=FYYYT(L) 
1330  NEXT  L 
1340  RETURN 
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ABSTRACT 
The  objective  of  this  work  was  the  evaluation  of  pseudospectral 
methods  for  their  computational  efficiency  and  applicablity  to  problems 
in  transport  phenomena.   This  was  accomplished  by  application  of 
Chebyshev  pseudospectral  methods  to  problems  from  each  major  area  of 
transport  phenomena.  The  evaluation  of  the  pseudospectral  method  was 
based  upon  comparison  to  the  analytic  solution,  if  available,  or  finite 
difference  approximation  otherwise.   It  was  found  that  discontinuities 
in  the  solution  domain  can  result  in  serious  deviations  from  the  correct 
solution;   for  example,  the  temperature  discontinuity  in  thermal 
entrance  length  problems  led  to  the  propagation  of  error  within  the 
solution.  Digital  filtering  was  used  successfully  to  damp  out 
oscillatory  behavior  in  all  cases  studied. 


